Dynamics and Control of Nonlinear Systems 


A Thesis Submitted 

in Partial Fulfillment of the Requirements 
for the Degree of 
Doctor of Philosophy 


by 

Rahul Konnur 


to the 

Department of Chemical Engineering 

Indian Institute of Technology, Kanpur 
April, 1997 



Certificate 



It is certified that the work contained in this thesis entitled Dynamics 
and Control of Nonlinear Systems by Rahul Konnur has been carried 
out under my supervision and that this work has not been submitted 
elsewhere for a degree. 



Associate Professor 
Department of Chemical Engineering 
Indian Institute of Technology 
Kanpur. 208 016 (India) 


April, 1997 



“9 JUL WVfl 



UBHUfiv 

• '• T - KAMPUfl 


M* A 


*25675 




Synopsis 


In this thesis we study the open-loop and closed-loop response of nonlinear sys- 
tems. In the first part we investigate the open-loop behavior of two nonlinear reacting 
systems using bifurcation theory and singularity theory. In the second part we apply 
results from differential-geometry based nonlinear control theory to analyse the prob- 
lem of sychronization of chaotic systems. A survey of the literature and the scope and 
motivation of the thesis are outlined in chapter 1. 

In the second chapter we discuss the dynamics of a discretely forced batch reactor 
system. Each cycle of operation of the reactor consists of the following steps. In the 
first step, the reactants are charged instantaneously into the reactor and the reaction is 
allowed to occur in batch mode for a fixed time P. At the end of this time, the reactor 
contents are emptied partially and the reactor volume is made up by the addition of 
fresh reactant. Both the discharge of the reactor contents and the recharging of the 
fresh feed are assumed to occur instantaneously. The above steps are repeated in the 
next cycle of operation. Hence the operation of the reactor is periodic, with period P. 
Such an operation emulates the batch process in the limit of P — > oo and the CSTR 
as P — )■ 0. 

We consider the reaction occurring to be the first order, exothermic and irreversible 
reaction A -4 B. The reactor dynamics in the two limiting cases of forced operation, 
viz., the batch reactor and the CSTR, are well understood. This investigation enables 
us to study the transition of the dynamics of this discretely forced reactor system from 
that of the batch reactor to the CSTR. We have used the shooting method and the 
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Fig. 1. Schematic diagram of the reactor-separator system with recycle. 


continuation method to compute the periodic solutions and generate the bifurcation 
diagrams of the system respectively. Stability of the periodic solutions has been com- 
puted using Floquet theory. The stroboscopic map and the Lyapunov exponents have 
been used to characterize the nature of the dynamic behavior of the system. Extensive 
s im ulations indicate that the system can exhibit a wide variety of complex periodic 
and aperiodic behavior. These include the period doubling, intermittency and quasi- 
periodic routes to chaos, crisis and mixed-mode oscillations. 

In chapter 3 we investigate the steady-state multiplicity behavior of a reactor- 
separator system with recycle (Fig. 1). We consider the reactor to be a CSTR again 
sustaining the first order, exothermic and irreversible reaction A -» B. The separator 
is considered to be a single stage flash. It separates the reactor effluent stream into a 
reactant rich liquid stream and a product rich vapor stream. The reactant rich liquid 
stream from the separator is recycled back to the CSTR, while the product rich vapor 
stream is withdrawn from the system. 

We study the effect of different sets of design specifications on the steady-state 
behavior of the system. The objective of this work is to investigate the steady-state 
multiplicity features and obtain a qualitative picture of the dynamic behavior that 
can be expected for different design specifications of this system. The system selected 
is such that it permits the application of singularity theory for the determination 
of the steady-state multiplicity features. For the coupled system, the assumption of 
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vapor-liquid equilibrium separation imposes constraints on the operation of the coupled 
system. These constraints influence the multiplicity features of the coupled system. 
Application of singularity theory enables us to obtain a comprehensive picture of the 
steady-state multiplicity features of the system in the presence of these constraints. 
The bifurcation diagram is used to represent the multiplicity features of the system. It 
is shown that for some sets of design specifications, there can exist regions in parameter 
space, where there are no steady-states of the coupled system. 

In the next two chapters we discuss the application of results from nonlinear control 
theory for the analysis of the phenomenon of synchronization in chaotic systems. Two 
systems are said to synchronize when they evolve along identical trajectories in phase 
space. 

Chaotic systems exhibit a sensitive dependence to initial conditions, i.e., two iden- 
tical chaotic systems starting from nearby initial conditions, evolve along exponentially 
diverging trajectories. This sensitive dependence on initial conditions associated with 
chaotic evolution means that two identical chaotic systems are not expected to syn- 
chronize with each other. However, it has been recently shown that it is possible to 
achieve synchronization in chaotic systems also. Several different techniques have been 
proposed for this purpose. 

In chapter 4 we present a unified framework for the analysis of the different tech- 
niques for synchronizing chaotic systems. We employ recent results from differential- 
geometry based nonlinear control theory for this purpose. It is shown that the different 
synchronization techniques are a different interpretation of perfect control of a sys- 
tem along a desired chaotic trajectory. By perfect control, we mean that the system 
(plant) output tracks its desired trajectory exactly for all t > 0. We illustrate the 
connection between synchronization and perfect control on the Rossler system. 

It is not possible to synchronize several systems using many of the techniques cur- 
rently available. The Rossler system is an example of such a system. These systems 
are classified as non-minimum phase systems in the control literature. These are the 
systems having zeros in the right-half plane. In the first part of the fifth chapter we em- 
ploy a control scheme based on bifurcation theory for achieving chaotic synchronization 
in the Rossler system. 

The synchronization problem is formulated as a problem of regulating the output 



of a system using state feedback. The reference output tracking signal is generated 
by a no nlin ear external system. This external system is selected such that the com- 
posite system consisting of the plant and the external system, has atleast one zero on 
the imaginary axis. This permits application of the center manifold theorem for the 
derivation of the control law for regulating the output of the composite system. 

Input multiplicity is a situation where different inputs can give the same output. 
A result of this is that the controlled system states cannot be stabilized at their de- 
sired steady-state. Systems exhibiting input multiplicity appear frequently in chemical 
engineering. A major limitation of many of the control techniques based on differential- 
geometry is their limited use for the control of systems exhibiting input multiplicity. In 
the second part of the fifth chapter we use the output regulation technique mentioned 
above to control an isothermal CSTR system, sustaining an autocatalytic reaction and 
exhibiting input multiplicity, around its steady-state. 
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Chapter 1 
Introduction 


In this dissertation we focus our attention on the behavior of nonlinear systems. We 
have chosen different model systems for our investigations and studied their open-loop 
and closed-loop characteristics. In. chapter 2 we discuss the dynamic behavior exhibited 
by a discretely forced batch reactor. Chapter 3 is concerned with the analysis of the 
steady-state characteristics of a coupled reactor-separator network. We discuss the 
effect of different design specifications on the behavior of the system. In these two 
chapters we concern ourselves with the open-loop response of the respective systems. 

In the fourth and fifth chapters we discuss the closed-loop performance of chaotic 
systems. In the fourth chapter we revisit the problem of designing a feedback controller 
with the objective of tracking a reference chaotic trajectory. We show that the study of 
this nonlinear control problem helps us to develop a unified framework for the analysis 
of the phenomenon of sychronization of chaotic systems. This technique is suitable 
only for a certain class of systems. In the fifth chapter we develop and apply a more 
general nonlinear control technique based on bifurcation theory for controlling a chaotic 
system. This technique is also employed for the control of a reactor system exhibiting 
input multiplicity. 

In the remainder of this chapter we discuss the literature related to each of the 
above problems mentioned above. We also provide the motivation for the investigations 
undertaken in each of the following chapters. 

The steady-state and dynamic behavior of the continuous stirred tank reactor 
(CSTR) has been the focus of several investigations over the last four decades. Even 
before the fifties, it was known that a few reaction systems could exhibit multiple 
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steady-states. Such behavior was first reported by Liljenroth in l01S[lj. Later Van 
Heerden[2] showed for the first time that the CSTR could exhibit similar phenomena. 
He determined the conditions for the existence of multiple stead\ -states in the i (Motor 
by considering the dependence of the rates of heat generation and heat removal on the 
reactor temperature. The subsequent interval has seen atleast, four bursts of activity 
aimed at understanding the behavior of this system. 

In the first phase, the emphasis was on developing a mathematical framework for 
understanding the uniqueness and multiplicity features of the CSTR. Most investiga- 
tions considered the exothermic, irreversible reaction scheme A -4 B in a non-adiabatic 
CSTR. Bilous and Amundson[3] applied Lyapunov theory to cam’ out a linear stability 
analysis for this system. Their investigations revealed the existence of different types 
of steady-states for this system. They applied the concept of phase plane diagrams for 
explaining the nature of approach to the steady-state. 

In the mid-seventies, the pioneering work of Poore and co-workers[4~6] resulted in 
a second burst of activity in investigations into the multiplicity features of the CSTR. 
This effort was mainly directed towards systematically characterizing the different pos- 
sible types of steady-state multiplicity and dynamic features of this system. Bifurca- 
tion theory was used for this purpose. They introduced the bifurcation diagram for 
graphically representing the dependence of a state variable on a independent control 
parameter. 

Multiplicity features in the context of non-isothermal reactions arise due to the 
presence of thermal feedback. This idea goes back to the work of Van Heerden[2]. 
Another class of reactions which have been widely investigated are the autocatalytic 
reactions. Here the multiplicity features arise because of feedback due to autocatalysis. 
Such autocatalytic mechanisms are of importance in catalytic and biochemical reaction 
systems. Gray and Scott [7, 8] have thoroughly investigated the multiplicity features 
of different types of autocatalytic reaction schemes in the CSTR. They showed that 
the autocatalytic mechanism by itself is sufficient for a system to exhibit a variety of 
multiple steady-state and dynamic behavior. 

In the early eighties, the focus of the investigations shifted towards understanding 
the dynamic behavior of the CSTR. This was an outcome of the discovery that several 
simple nonlinear systems could exhibit complex dynamic behavior. Jorgensen and 
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Aris[9] demonstrated the richness of the dynamic behavior of the reaction A -> B — » C 
in a CSTR. Such complex dynamic behavior has also been observed in non-isothermal 
autocatalytic reactions in a CSTR[10]. 

Some studies on the complex dynamics in CSTR’s have considered the reaction 
A — y B in a periodically forced system. It has been shown that periodic variation of the 
coolant temperature or the reactant flow rate can cause the system to exhibit complex 
dynamic behavior. Sincic and Bailey [11] have studied the dynamics of the first order, 
exothermic reaction in a periodically forced reactor. They investigated the dynamics 
under very fast forcing, very slow forcing and at intermediate values of the forcing 
frequency. They showed that it was possible for the system to exhibit sub-harmonic 
and quasi-periodic oscillations at intermediate ranges of the forcing frequency. 

Several studies have been carried out to explain the mechanism through which such 
complex behavior can arise in the context of the CSTR. Aris and coworkers[12-15] have 
investigated the bifurcation features of different reaction schemes. Their objective was 
to study the effect of the bifurcation and multiplicity features of the unforced system 
on the system with forcing. They suggested mechanisms for explaining the occurrence 
of complex dynamic behavior in the periodically forced CSTR. Studies on verifying the 
universal features of such complex dynamics have also been reported in the literature. 
For example, Mankin and Hudson[16] have demonstrated the existence of a period 
doubling route to chaos and bistability for this system. 

Other than its importance from the point of view of nonlinear dynamics, periodic 
operation of chemical reactors has been widely studied in the literature. In some cases, 
it may be advantageous to operate the reactor in a periodic manner. For example, 
under certain conditions of operation, the reactor may have multiple steady-states. 
Some of these states may possess undesirable features. For instance, a steady state 
could have a low conversion, and another could be at a very high temperature, for 
which a reactor may not be designed. Alternatively, a steady-state could be unstable. 
In such cases, periodic operation of the system can allow one to overcome many of 
the drawbacks associated with the presence of multiple steady-states while combining 
their advantages. It has been reported that use of periodically forced reactors can 
result in higher reaction rates, improved selectivity and better controllability charac- 
teristics^?]. The parameters which can be varied periodically include the feed flow 
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rate, feed concentration and the coolant temperature. 

The above studies are based on the continuous forcing of a system. Another mode 
of a periodically forced operation of a reactor is discrete forcing, i.e., forcing at discrete 
times. Codell and Engel[18] were the first to investigate the performance of such a 
system. They considered the isothermal and adiabatic operations of a control led cycled 
tank reactor (CCTR). The CCTR is a batch reactor in which a pre-determined fraction 
of the reactor contents are withdrawn and replaced with fresh feed after fixed intervals 
of time. They showed that for some reactions, the performance of the CCTR could be 
better than a CSTR and a plug flow reactor. Ausikaitis and Engel [19] experimentally 
studied the reaction between sodium thiosulphate and hydrogen peroxide in a CCTR. 
They showed that the CCTR offers better controllability characteristics in the case 
of exothermic reactions than a CSTR. Lin and Wu[20] have carried out a modelling 
and experimental study of the reaction between sodium thiosulphate and hydrogen 
peroxide in an adiabatic CCTR. They have considered a more general situation of 
the case considered by Ausikaitis and Engel[19], with the objective of predicting the 
approach to the steady periodic state. Another of their objectives was to compare the 
performance in the CCTR with that in the CSTR and the plug-flow reactor. 

Kubickova et a/. [21] have investigated the nonlinear dynamics of the first order, 
exothermic reaction in a non-adiabatic fed-batch reactor (FBR). This mode of opera- 
tion of a reactor is an idealisation of the system considered by Codell and Engel [18], 
as they neglect the time of addition and withdrawal of fresh reactants and products 
respectively. They employed the arc-length continuation method to obtain the depen- 
dence of the periodic behavior of the system on the residence t im e in the reactor. They 
showed that the system could exhibit a variety of complex dynamic behavior. 

In the next chapter, we investigate a more general model of the system considered 
by Kubickova et al. [21]. Operation of the FBR is considered to be to a mode between 
the batch and continuous modes of operation. The dynamics of the batch reactor are 
uninteresting as the reactor always proceeds towards a state of complete conversion, 
as we consider the reaction to be irreversible. The CSTR, on the contrary, can exhibit 
different types of steady-state and dynamic behaviour, i.e., hysteresis, isolas and limit 
cycles. For this reaction, the dynamics of the CSTR are well known. Our objective is 
to study the dynamics of the FBR, and, to compare it with that of a CSTR operating 
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under the same conditions. This would enable us to see how the dynamics of the FBR 
evolves from that of the CSTR. In particular, we study how the different bifurcations 
and multiplicity features of the CSTR manifest themselves in the FBR. With this 
in mind, we have used the arc-length continuation method to obtain the bifurcation 
diagrams of the FBR. The stability of the periodic solutions have been determined 
using Floquet theory. Extensive simulations have been carried out to corroborate the 
predictions of the stability theory. 

The research on understanding the behavior of chemical reactors has been aided by 
the application of novel techniques from mathematics. The final burst of activity on 
investigations into the steady-state multiplicity features of chemical reactors followed 
the work of Balakotaiah and Luss [22-24]. In a series of papers, they demonstrated 
the application and usefulness of singularity theory in analysing the steady-state mul- 
tiplicity features of several lumped parameter systems. Since then singularity theory 
has been shown to be a powerful tool specially for the steady-state analysis of systems 
containing a large number of parameters. Application of singularity theory has made it 
possible to obtain the steady-state multiplicity features of several single reactor systems 
with different types of nonlinearities. 

Single units seldom occur in the process industries. It is typical to have several 
different types of units linked together by heat or mass exchange. Introduction of 
coupling through material and heat recycle in chemical processes alters the dynamics 
from that of single units. Therefore, eventhough the behavior of single units is well 
known, it is very difficult to predict the behavior of the coupled system. General 
features of the dynamic behavior of integrated plants have recently been reviewed by 
Morud and Skogestad[25, 26]. They have also have pointed out the need of obtaining 
some ’rules of thumb’ which could be used as indicators for classifying different types 
of behavior of the coupled system. 

Recent studies have analysed the steady-state and dynamic behavior of coupled 
reactor systems[10, 27]. These studies show that even when a model of the system 
is available, a comprehensive analysis of all the possible types of behavior of such, 
systems is very difficult. In principle, currently available mathematical tools such as 
bifurcation theory and singularity theory, can be used to analyse a plant. However, 
these tools cannot be used to obtain a complete picture of the behavior of anything 



but the simplest of such systems. 

Luyben and co-workers[28-30] have studied the dynamics and control of binary and 
ternary recycle systems. They consider a process consisting of an isothermal CSTR 
coupled through recycle to a distillation column. The effect of operating the system 
under different design specifications on the steady-state economics and the controllabil- 
ity of the process were discussed. Starting with a given set of design parameters, their 
approach consists of developing a design procedure for the operation of the distillation 
column. The effect of changing the design parameters on the system variables was then 
studied. Based on their results, they have conjectured about the advantages of working 
with certain sets of design parameters. However, the. complexity of the system in terms 
of the large number of variables, makes it impossible to obtain any firm results. In a 
later paper[31], Luyben analysed the effect of different control configurations on the 
snowball effect. The snowball effect is associated with a very large change in the recycle 
flow rate for a small change in the load variable. However, they did not provide any 
explanation for the origin of this phenomenon and ways for preventing its occurrence. 
In a later series of papers, they have extended their investigations to cover different 
cases of ternary systems [3 2-34]. In these studies, they consider several model reaction 
schemes, multiple separator units and multiple recycle streams. 

In the third chapter of this dissertation, we study the steady-state behavior of a 
reactor-separator system with recycle. As mentioned earlier, the overall behavior of 
systems with material or heat recycle is very different from the behavior of the indi- 
vidual units. Further, it has been claimed that plant interconnections may introduce 
fundamental limitations in the performance of any control system[25, 26]. Knowledge 
of the performance of such plants is important for control design. Hence it is impor- 
tant to have a thorough knowledge of the steady-state and dynamic behavior of recycle 
systems. With this in mind, we have investigated the steady-state multiplicity features 
of a model reactor-separator system with recycle. The system selected is such that it 
permits the application of singularity theory [35]. Application of this technique enables 
us to obtain a comprehensive picture of the steady-state multiplicity features of this 
system. It is shown that this system exhibits some unusual features. For example, 
we find that the system cannot have any steady-states for some sets of the oeprating 
parameters. We provide a physical explanation for the origin of these unusual behav- 



7 


ior. The objective of this work is to obtain a qualitative picture about the types of 
steady-state behavior that can be expected for different sets of specifications of this 
system. For a plant with recycle streams, this study provides a staxting point for a 
more detailed analysis and for understanding the different physical interactions. 

Over the past decade, there have been major developments in the applications of 
nonlinear control schemes based on differential-geometry for the control of chemical 
engineering systems[36]. Extensive experimental and simulation studies on many sys- 
tems, using several control techniques have been reported in the literature. In chapters 
4 and 5 of this thesis, we show the applicability of results from nonlinear control theory 
to a problem which has been attracting widespread attention in the physics literature. 

Since 1990, synchronization of chaos has been a topic of immense interest and 
puzzlement in the physics community. Two systems are said to synchronize if, after 
starting from different initial conditions, they evolve along the same trajectory in phase 
space. Synchronization between two identical chaotic systems has been considered to 
be an unlikely goal, because chaos is characterized by a sensitive dependence on ini- 
tial conditions. However, Pecora and Carroll[37, 38] showed that two chaotic systems 
could be synchronized with each other. Within the last two years, two different types 
of synchronization have been reported in the literature. Rulkov et al. [39] and Kocarev 
and Parlitz[40] have studied the type of synchronization which they refer to as gen- 
eralized synchronization. This type of synchronization is characterized by a nonlinear 
relationship between the drive and response system variables. Kocarev and Parlitz and 
their co-workers [41, 42] have reported a powerful new technique for achieving synchro- 
nization in chaotic systems! Their procedure, an active-passive decomposition (APD) 
technique involves creating a new system by defining a new variable as a function of 
the system variables and using this to drive the modified response system. 

The importance of synchronization stems from several potential applications [40]. 
For example, it has been shown that synchronization can be used in secret communi- 
cations and cryptography. Another potential application of synchronization is in the 
understanding of neural processes, where such a mechanism can be used to simulate 
neural behavior. There has also been much speculation in the literature about the role 
of synchronization governing the smooth functioning of the heart. 

In the fourth chapter, we present an approach for the analysis of the different types 
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of chaotic synchronization in a unified framework. This approach results from t he use of 

results from nonlinear&ntrol theory, which have been employed by chemical engineers 
in recent years. ' 3 


Systems exhibiting input multiplicity are frequently encountered in chemical en- 
gineering. Input-multiplicity is a situation where different inputs can give the same 
output. Balakotaiah and Luss[45] have applied singularity theory to classify a param- 
eter space into regions where it is possible to encounter input multiplicity. Koppel[43 
44] has carried out a detailed analysis of the effects of input multiplicity on different 
model control systems. However there has not been much progress in the development 
o control techmques for use on such systems. This is especially true for several of the 
me hods based on differential-geometry. A limitation of many of these control tech- 
niques us that the, can be used only for the control of minimum-phase systems. Hence 
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Chapter 2 

Nonlinear Dynamics of a 
Fed-Batch Reactor (FBR) 
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Period of fed-batch operation. 
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2.1 Introduction 

In the previous chapter we reviewed the literature on the nonlinear dynamics of chem- 
ical reactors. Most of these studies have concentrated on investigating the dynamic 
behavior e xhi bited by the CSTR. The dynamics of several reaction schemes in (i) the 
autonomous CSTR and (ii) the periodically forced CSTR have been the subject of 
several detailed studies[10]. In comparison, the nonlinear dynamics of discretely forced 
systems have not been widely studied. 

Discretely forced systems occur in a variety of contexts. In the area of medical 
sciences, periodic variation of drug concentration in the body is an example of such a 
system. Intake of the drug results in a sudden increase in its concentration in the body. 
Subsequently, ingestion of the drug results in a gradual decrease in its concentration. A 
similar situation also prevails in the periodic addition of fertilizers in farming. Discrete 
forcing is also widely used in the manufacture of various products involving biochemical 
reactions. Here, it is often desirable to add the substrate to a batch bio-reactor as the 
reaction proceeds. This addition is often done in discrete pulses, in order to optimize 
the reactor performance. 

In this chapter we investigate the dynamic behavior of the first-order, irreversible 
and exothermic reaction in a non-adiabatic fed-batch reactor (FBR). The FBR is an 
example of a discretely forced batch system. The mode of operation of the FBR is in 
between the batch mode and the continuous mode. The dynamics of the batch reactor 
are uninteresting, as the reaction always proceeds in the direction of increasing con- 
version. The CSTR, on the other hand, exhibits a rich variety of dynamic behavior[6]. 
We wish to compare the dynamics of the FBR with that of a CSTR operating under 
the same conditions. In particular, we will be interested in investigating how the mul- 
tiplicity, stability features and the nature of approach to the steady-state in the CSTR, 
influence the dynamics of the FBR. This would enable us to see how the dynamics of 
the FBR evolve from that of the corresponding CSTR. Therefore, our broad objective 
will be to investigate how the multiplicity features of the FBR depend on one of the 
system parameters. For this purpose, we have used an arc-length continuation method 
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to compute the bifurcation diagrams of the FBR. The shooting method has been used 
to compute the periodic solutions of the FBR model. The stability of the solutions 
have been determined using Floquet theory. Extensive simulations have been carried 
out to bring out the different types of complex dynamic behavior exhibited by the 
FBR. 

We start by describing the FBR model. This is followed by a brief description of 
the method of solution. The results axe presented in the third section and we conclude 
with a discussion of the results. 

2.2 Description of the Model 

The FBR discussed in this chapter is a periodically forced batch reactor. The operation 
consists of the following steps [21]: 

1. The reactor is filled up to the desired volume with the reactants at the start of a 
operating cycle, 

2. The reactor is well stirred and the reaction is carried out in the batch mode for 
time P (reaction step), 

3. At the end of the reaction time, a fixed volume fraction R of the reactor volume 
is withdrawn over a time interval t d (discharge step) , 

4. The reactor volume is made up to the final volume by the addition of fresh feed 
over a time t c (recharge step), 

5. The steps 2, 3 and 4 are repeated. 

The FBR described above is different from that widely used for carrying out bio- 
chemical reactions. In the biochemical context, a FBR is a semi-batch reactor. 

Throughout this study, we assume for simplicity that t c — t d — 0. This implies that 
steps 3 and 4 are instantaneous. The system we consider is the first order, exothermic, 
irreversible reaction A B in a non-adiabatic FBR. Thus our model is the most 
general possible for this model reaction. The evolution of the concentration and the 
temperature in the reactor occurs during step 2. This evolution is modelled by the 
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following set 
respectively 


of dimensionless equations, representing the material and energy balances 


§ = k 0 r(l-C)e T = /, 


( 2 . 1 ) 


§ = Bk 0 r(l - C)e T - I3 0 tT = h 
with the initial condition: t = 0: C = 0, T = 0. 

Here we have assumed that the concentration and temperature characterising the 
initial state of the reactor (at the beginning of the reaction step of the first cycle) are 
equal to that of the fresh feed used for recharging the reactor. In the above equations, 
B, P 0 and k 0 represent the dimensionless heat of reaction, the heat removal parameter 
and the reaction rate constant respectively. The parameter r (= P/R) is representative 
of the residence time in the reactor. We have used r as the characteristic time scale to 
define the dimensionless time t. For simplicity, we have made the positive exponential 
approximation in the derivation of the above equations. 

Equation (2.1) is valid only during the interval of batch operation. At the end of 
operation in the batch mode (of interval P actual time units), the discharge of products 
and the recharge of fresh feed updates the values of C and T in the reactor. In terms 
of the dimensionless time t, this is equivalent to updating after every integral multiple 
of R dimensionless time units. Therefore, at the end of each cycle of batch operation 
(of R dimensionless time units), the discharge of products and the recharge of fresh 
feed updates the values of C and T in the reactor according to 


C t=R+ = (1 - R)C t=R - 


( 2 . 2 ) 


T t=R+ = (1 -R)T t=R . 

In the above equations, the subscripts — and + denote the values of the states before 
and after the updating respectively. Therefore, at the end of operation in the batch 
mode, there is a discontinuous variation of the conversion and temperature in the 
reactor governed by (2.2). 

The model of the FBR appears to be an autonomous system as the independent 
variable t does not explicitly occur in the model. This system is a periodically forced 
system and is forced with a period P . The evolution of the system is therefore, deter- 
mined not only by the current state of the system, but also by the time at which the 
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state is reached. Hence, if a state is reached at the end of the reaction step, it changes 
discontinuously according to (2.2), else its evolution is governed by (2.1). This is a 
characteristic of a non-autonomous system. 

The FBR as described above, tends to a batch process in the limit of R — ¥ 0. In 
this limit the reaction approaches a state of complete conversion with increasing time. 
This follows since the reaction is irreversible. For R = 1, we are at the limit of batch 
operation which is repeated every P time units. These two limits are uninteresting as 
far as the dynamic behavior of the reactor is concerned. However, for 0 < R < 1, we 
have an open system which exchanges mass and heat with the environment periodically. 
We therefore, analyze the behavior of the system for R € (0, 1). For finite P and 
R € (0,1), the mean residence time r in the reactor can be defined as P/R. The 
dynamic characteristics of the FBR can hence be compared with that of a CSTR of 
the same residence time, for identical values of other system parameters. The true 
limit of the CSTR is reached in the limit P -> 0, R — » 0, such that 

p 

lim — = t ( finite ) 
r-+o,p-+o R w ’ 

The CSTR of residence time r c , which is equivalent to the FBR described above can 
be represented by the following set of equations[6] 

f = -C + k 0 r c (l - C)e T 
f = -T + Bk 0 r c {l - C)e T - /3 0 r c T 

2.3 Method of Solution 

The basic state of the CSTR is a steady-state. In the FBR, the equivalent state is a 
periodic state of the reactor which has a period P. Hence this is referred to as a IP 
state. A state which is periodic and has a period mP with m > 1, is denoted as an 
mP state. This corresponds to sub-harmonic oscillations of the FBR. As mentioned 
earlier, the mP states (jn P 1) are actually obtained as mR states, since we shall be 
using the dimensionless set of equations (2.1) and (2.2). However, we employ the mP 
notation to identify these periodic states of the FBR. 
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In the bifurcation study of a system, one investigates how a dependent, variable 
varies with a control parameter. This is called the bifurcation paiaxnetei . hui example, 
in the classical approach employed for understanding the bifurcation characteristics of 
the CSTR, one studies the dependence of the steady-state conversion or temperature 
on the bifurcation parameter (the residence time) in the reactor[6j. In this work we 
determine the multiplicity and stability features of the IP state of the FBR. The 
par am eter P is considered to be the bifurcation parameter. Changing P is equivalent 
to changing the residence time in the FBR. 

We begin by describing the shooting method which has been used to generate the 
IP solutions of (2.1) subject to (2.2). Later, we discuss the arc-length continuation 
method which has been used to generate the dependence of C and T on P. This is 
followed by a discussion on the determination of stability of the IP solution. Details 
of all these techniques axe available in [48]. 


2.3.1 Shooting Method 

The IP solution of (2.1) subject to (2.2) can be determined by integrating (2.1). once 
the initial conditions C 0 and T 0 of the periodic state have been determined. It is 
possible to determine these unknowns by formulating and solving a boundary-value 
problem using the shooting method. For a IP solution, we require C I{ f = Co and 
T r+ = T 0 . This condition is represented as 


G 1 (C 0 ,C r+ ,T o ) = C r+ -C o = 0 


G 2 (C q , T r+ , T 0 ) = T r+ - T 0 = 0 



The equations (2.3) are a set of algebraic equations and these can be solved for the 
unknowns Co and To using the Newton-Raphson method. The values of the variables 
before and after updating at t = R (denoted by the subscripts — and 4- respectively) 
are dependent on the values at the beginning of the cycle (denoted by a subscript 0). 
These are determined by integrating (2.1) with the initial conditions (C 0 ,T 0 ). This 
is followed by imposing the updating condition (2.2). The (n + 1)* iterate of the 
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unknowns is computed from the n th iterate using 



n+1 


n 

dCji+ -| 8Cr+ 

n- 1 
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dC 0 1 dT 0 


’Co' 

. To _ 
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dC 0 dT 0 J 




(2.4) 


The variations in C R+ , T R+ with Co and T 0 are obtained by solving the following set 
of four linear ordinary differential equations 


A (9C_\ _ (9jA ( dC_ \ , (djA ( er \ 

dt\dCoJ \dC J \dCoJ ^ \dT J \dCoJ 

A (dc\ _ (ejA (ac\ , (ajA far) 
dt \dTo) \9C J \dToJ ^ \&T ) [dToJ 

A(AL) = (Ml) (AC) , /sM (jr.\ 
dt \dC 0 J V. dC ) V dCo J ' \ dT ) \ dCo ) 

A(ML) _ (Ml) (8C\ + (Mi) (AT) 
dt\dT 0 ) \dC J \dToJ ^ \dT J \dToJ 


with the initial conditions 


At t = 0, 


ac_ _ i 9C-—n ML _ n ML — 

dCo ~ ’ dTo ~~ u ’ dCo ~ > dT 0 ~ 


(2.5) 


and 

—kore T k Q T(l — C)e T 

-Bk 0 re T Bk 0 r(l - C)e T - for _ 

The system of six equations (2.1) and (2.5) are integrated simultaneously. Depen- 
dence of these variables on the updating condition is accounted for, by imposing the 
updating condition on these variables. The values of the derivatives in the matrix 
of (2.4) are the values of the dependent variables in (2.5), at t = R+. We can thus 
obtain the solution C 0 , T 0 using an iterative procedure (using (2.4)) for a fixed value 
of the system parameters B, k 0 , P, R and fy. Once C 0 , T 0 have been obtained, the 
IP solution is obtained by integrating (2.1) using. (2.2). Equation (2.3) is a periodic 
boundary condition. This ensures that the dependent variables C and T are periodic 


Ml 

Ml 1 

dC 

d T 

Mi 

Mi 

dC 

dT J 
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with a period R. For the cases where multiple IP solutions exist for (2. 1 ), the solution 
to which we converge upon, depends on the initial guess ot Cq and 

Having determined the solution for a fixed set of parameters, we fix all the param- 
eters except P. We study how the IP solution changes as we vary /’ - the bifurcation 
parameter. The arc-length continuation method has been used to obtain the depen- 
dence of the IP solutions of the system (2.1) on P. This method is described next. 

2.3.2 Arc-length Continuation Method 

It is well known that for the first order reaction, the bifurcation diagram of the (,'STR, 
can consist of different branches. Hence we may expect that for the FBR, the depen- 
dence of the IP state on the bifurcation parameter P can consist of various branches. A 
technique which can be used to determine the dependence of the dependent, variable on 
an independent variable is the arc-length continuation method[48]. We have employed 
this technique for generating the dependence of the IP states of the FBR on P in a 
smooth manner. The method involves introducing a new parameter, tin; arc-length of 
the solution curve denoted by s, such that the system has a smooth set of solutions in 
(C, T, P) space. The method is based on starting with a solution for a given value of 
P, say P 0 . We use this solution as an initial guess to determine the solution at a new 
value P + SP in the neighborhood of the initial point P 0 . We now describe the details 
of this technique. 

Let s denote the arc-length along the solution curve. Having determined the IP 
solution of (2.1) subject to (2.2) at an initial point P 0 , we denote that point by s = 
0. Instead of studying the dependence of C(t) and T(t) on P, the problem is now 
transformed to that of studying the dependence of C(t), T(t) and P on the arc-length 
s, i.e., we are interested in determining the evolution of C{t), T{t) and P as $ varies. 
The IP solution of (2.1) subject to (2.2) is uniquely determined by its initial conditions 
C 0 and T 0 , for a given P. Therefore, the arc-length continuation procedure reduces to 
a problem of determining how C 0 , T 0 and P change, as we proceed along the arc-length 
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s. Differentiating (2.3) with respect to s, we have 


* = (») («0 + (fit) (f 1 ) + (») (f ) = o 

( 2 . 6 ) 

* = (») w + (») (f 1 ) + m (f ) = o 

The coefficients and occurring in the derivatives of Gj and G 2 in (2.6) are 
obtained by integrating the following equations 


ft (SO - (so) (S) + (It) (ff) + (It) 

(2.7) 

1(f) = (»)(») + (#)( ») + (») 


with the initial conditions: at t = 0: §£ = 0 and — = 0, followed by imposing the 
updating condition. 

Equation (2.6) is a system of two linear homogenous equations in the three un- 
knowns ~[f- and A third equation is obtained using the normalisation condition 


from calculus 



From (2.6), we seek a solution of the form 


( 2 . 8 ) 


i = l. 2, k- 1, * + l, •••, n+ 1 (2.9) 

as as 

In the above equation, Xi = Co, x 2 = T 0 , x 3 = P and n = 2. We can thus determine 
fa and @2 from (2.6) and solve (2.8) for ^ to evaluate 


dxk 1 

is ±p+itfTitf 


( 2 . 10 ) 


The sign of the derivative ^ is given by the orientation of s along the solution branch. 
Choosing the variable x k and the sign of the solution to the above equation appropri- 
ately, enables us to traverse through the bifurcation points in a smooth manner. 

We integrate (2.1), (2.5) and (2.7) simultaneously, with the updating conditions 
(2.2). Equation (2.9) is then solved for fa and fa using Gaussian elimination. We then 
obtain estimates of g, *§■ and from (2.9) and (2.10). This system of equations 
with the initial condition 
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At s = 0: P = P 0 , C = Co and T = T„ 

is integrated using the Euler method. This is the predictor step which yields new 
approximate values of P, Co, Tq. These values of Cq and Fq are then used as initial 
guesses to solve (2.1) and (2.3) with the new value of P , using the shooting method 
described earlier. This shooting method step is the corrector step, used to obtain 
the accurate IP solution with the new value of P. The Gaussian elimination. Euler 
integration and the shooting method are repeated with the new value of P computed 
above. This procedure eventually enables us to determine the dependence of the IP 
solution on P in a continuous manner. 


2.3.3 Stability 


IP solutions of (2.1) subject to (2.2) computed using the shooting method can be stable 
or unstable. The solution on which we converge is solely determined by the initial 
guesses of Co and T 0 used in the scheme. Floquet theory[48] is used to determine 
the stability of the IP solution. After a IP solution is determined, its stability is 
determined by calculating the magnitudes of the eigenvalues of the variational matrix 


' 1 
dCo dTo 


9Tr± 3Tr± 

L dC 0 dT 0 J 


( 2 . 11 ) 


The solution is stable if both the eigenvalues of the above matrix lie inside the unit 
circle in the complex plane. As the bifurcation parameter P is varied, the eigenvalues 
can cross the unit circle in three ways. Therefore, there are three possible ways in 
which the IP state can be de-stabilized. These are 


1 . One of the eigenvalues leaves the unit circle through — 1 . This corresponds to 
a period doubling bifurcation. The IP state of the FBR becomes unstable and 
another state with twice the original period branches off at this point. 

2. An eigenvalue leaves the unit circle through +1. In the case of the FBR, this 
instability corresponds to the saddle-node bifurcation. At this bifurcation point, 
a IP state changes its stability and another IP state emanates from this point. 
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3. A complex-conjugate pair of eigenvalues leaves the unit circle. This corresponds 
to a torus bifurcation. Following this bifurcation, there is a change in the dimen- 
sion of the attractor. We can now visualise the system as evolving on the surface 
of a torus. 

2.3.4 Simulations 

Extensive simulations were carried out to obtain a complete picture of the dyn ami c 
behaviour of the fed-batch reactor. Gear’s method (Routine D02EBF of NAG library) 
was used for integrating (2.1). Integration was usually carried out till 5000 cycles of 
operation of the fed-batch reactor. 

2.3.5 Stroboscopic Map 

In general it is very difficult to identify the nature of the dynamics using time series 
data obtained from the simulations. In such cases, the stroboscopic map is useful for 
identifying the nature of evolution of the system. This is a standard technique used 
in the analysis of forced systems. This technique involves measuring one of the system 
states after definite intervals of time. The strobing is generally done at a frequency 
equal to the forcing frequency. Construction of the stroboscopic map therefore, involves 
plotting the value of a state variable X(t) (either C or T) versus X(t+R). In this work 
we have used the value of T at the end of reaction step, to generate the stroboscopic 
maps. The stroboscopic map for the FBR showing mP dynamics consists of m points. 
For a system undergoing quasi-periodic dynamics, the stroboscopic map yields a closed 
curve, indicating that the system evolves on the surface of a torus. A stroboscopic map 
in which the points are scattered is representative of chaotic dynamics. These points 
usually lie on dense segments, which is typical of deterministic chaos. 

2.4 Results and Discussion 

The dynamics of our model reaction in the CSTR has been thoroughly investigated by 
Uppal et al.[6]. They catalogued the different steady-state multiplicity features and the 
qualitative aspects of the dynamic behavior of the CSTR by classifying the parameter 
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Set 

B 

k o 

A) 

1 

14.0 

0.162 

3.0 

2 

8.0 

0.136 

1.0 

3 

6.88 

0.136 

0.72 


Table 2.1: Parameter values. 



Figure 2.1: Bifurcation diagram for the CSTR for parameter sets 1 and 2 of Table 2.1, 
( — ): stable steady-state; (•••): unstable steady-state; (x): Hopf bifurcation point. 

space into different regions. For the purpose of comparing the dynamic characteristics 
of the FBR with the CSTR, we have chosen three sets of parameter values for our 
investigation. These are listed in Table 2.1. 

The parameter values in set 1 are identical to those chosen by Kubickova et al. [21]. 
For this set, the CSTR has multiple steady-states for 0.00183 <r c < 0.5458. The upper 
steady-state branch is unstable in the region 0.5963 < r c < 0.9163. This instability 
arises because of a Hopf bifurcation and the system exhibits oscillatory behavior in this 
region. The dependence of the steady-state conversion of the CSTR on r c is shown in 
Fig. 2.1. Kubickova et oh[21] have shown that the FBR exhibits periodic oscillations 
of different periods for parameter values of this set under the restriction P = R. In 
the present study, we generate the bifurcation diagram of the IP state of the FBR 
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Figure 2.2: Bifurcation diagram of the IP state for the FBR, parameter set 1, 
R = 0.573875, ( — ): stable IP state; (•••) : unstable IP state; (x): Period doubling 
bifurcation; (•): Torus bifurcation point. 

by varying P, for a fixed value of R and the other parameters. This is equivalent to 
determining the effect of residence time on the IP state of the FBR. 

We concentrate only on the basic periodic behavior of the FBR. This is the periodic 
behavior with a period equal to the forcing period P. As mentioned earlier, we refer 
to this as a IP state. This IP state is analogous to the steady-state of the CSTR. We 
have shown the dependence of the IP state on P for R = 0.5738750 in Fig. 2.2. The 
ordinate in the figure is the conversion at the end of a batch cycle before the updating. 
This corresponds to the composition of the product withdrawn from the reactor. The 
figure is qualitatively similar to the bifurcation diagram of the CSTR (Fig. 2.1). There 
is a region of r where multiple IP states exist and a region where the upper IP state 
branch is unstable. The different features of Fig. 2.2 are discussed next in detail. 

For R = 0.573875, as we increase P from zero, the conversion increases till P = 
0.410027. At this point an eigenvalue of the variational matrix (2.11) moves out of the 
unit circle through +1 and the IP state loses its stability via a saddle-node bifurcation. 
This corresponds to the ignition point of the system. The IP state emerging as a result 
of this bifurcation is a saddle, as one eigenvalue of (2.11) lies inside the unit circle and 
the other is outside the unit circle in the complex plane. This unstable state exists 
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for P < 0.410027. This IP state exists till P = 0.00251, where another saddle-node 
bifurcation occurs. Here an eigenvalue of (2.11) re-enters the unit circle through +1. 
The IP state now regains its stability and it turns around. This point corresponds to 
the extinction point. The stable branch exists for P > 0.00251. As we continue along 
P on the top branch, an eigenvalue of (2.11) goes out of the unit circle through —1 at 
P = 0.3186. At this value of P, the IP state loses its stability in a period doubling 
bifurcation. As P is increased past this value, the stable state of the FBR is a periodic 
state with a period of twice the forcing period. We now continue along this stable 2P 
state and study successive bifurcations. As we increase P, we encounter a sequence 
of period doubling bifurcations. The first five period-doubling bifurcations points are 
given by the first five values of P in Table 2.2. At each of these values of P, an eigenvalue 
of the stable 2 m P state (m = 1 to m = 5),. existing to the left of the bifurcation point, 
moves out of the unit circle through —1. The distance between successive values of P 
at which period doubling bifurcations occur decreases very rapidly. This is consistent 
with the classical period doubling route to chaos[49]. This cascade of period doubling 
bifurcations continue till the bifurcation points accumulate at a value of P close to 
0.33513. We were unable to continue beyond the 32P state. This is because the 
bifurcations occur very rapidly. 

As we vary the parameter P from 0.33515 to 0.351, i.e., just beyond the regions 
of period doubling bifurcations, the FBR exhibits interesting dynamic behavior. The 
stroboscopic map for P = 0.33515 is shown in Fig. 2.3(a). The last 2000 points out 
of 5000 cycles of the simulation were used to generate the map. This stroboscopic 
portrait can be described as two small islets facing each other. Further, the strobed 
points oscillate in the neighbourhood of T = 4.5 and T = 6.6. This dynamic behavior is 
typical of the emergence of chaotic behavior from the classical period doubling route to 
chaos. The stroboscopic maps for P = 0.342 and P = 0.349 are shown in Figs. 2.3(b) 
and 2.3(c). This type of the stroboscopic portrait is typical of a system exhibiting 
chaotic dynamics. The growth and the subsequent merging of the islets as we increase 
P is clear from Figs. 2.3(a)-2.3(c). The three cases considered above show that the 
system exhibits chaotic behavior, for values of P just beyond the accumulation point of 
the period doubling cascade. This behavior ends suddenly for values of P just beyond 
0.351. For 0.351 < P < 0.410027, i.e., till the ignition point, the only stable state of 
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Nature of FBR dynamics 

P 

1P-2P 

0.3186 

2P-4P 

0.332955 

4P-8P 

0.334596 

8P-16P 

0.3350285 

16P-32P 

0.3351235 

Chaos 

0.33513-0.351 

Crisis (lower IP) 

0.351-0.410027 

Intermittency 

0.410027-0.411 

Mixed-mode oscillations (Table 2.3) 

0.415-0.70 

Torus bifurcation 

0.8214 


Table 2.2: Bifurcation points of the FBR and summary of the dynamic behavior for 
parameter set 1 of Table 2.1, R = 0.573875. 

the FBR is the IP state corresponding to the lower stable IP branch in Fig. 2.2. 

In this region the chaotic attractor loses its basin of attraction, i.e., points in phase 
space are no longer attracted to it. This is induced by the collision of the chaotic 
attractor with the upper unstable IP state. This phenomenon is termed as a boundary 
crisis [50, 51]. This can be seen from the stroboscopic map for P = 0.349 shown in 
Fig. 2.3(c), which shows the chaotic attractor about to collide with the fixed point. 
Destruction of the chaotic attractor results from its collision with the upper unstable 
IP state near P = 0.351. For values of P > 0.351, the system exhibits transient chaos, 
i.e., the trajectory stays in the vicinity of the chaotic attractor (which has disappeared 
now) for some time before being attracted to the lower IP state. The crisis point is a 
bifurcation point as beyond this point the chaotic attractor disappears. This bifurcation 
involves two non-local attractors and hence is an example of a global bifurcation. 

There is another sudden change in the dynamic behavior as we cross the value 
of P at the ignition point (P = 0.410027). For values of P > 0.410027, the lower 
IP state does not exist, and the upper IP state is unstable. The steady behavior of 
the system is such that it consists of a large number of small amplitude amplitude 
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T(t) 


Figure 2.3: Stroboscopic maps of the chaotic attractor showing the growth and collision 
of the chaotic attractor with the unstable fixed point (x), parameter set 1, R = 
0.573875, (a) P = 0.33515; (b) P = 0.342 and (c) P = 0.349. 
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oscillations which are occasionally interrupted by a few large amplitude oscillations. 
This behavior is characteristic of a system exhibiting Type I intermittency[52, 53], 
where we have bursts of noisy behavior interrupting an otherwise quiescent evolution 
of the system. Figure 2.4(a) shows the variation of T with dimensionless time for 
P = 0.4105. This is actually a 78P state characterized by a large number of small 
amplitudes peaks and a few big amplitude peaks. Figure 2.4(b) shows the variation 
of T with time for P = 0.4110. The stroboscopic maps of the attractors for these two 
values of P are shown in Figs. 2.4(c) and 2.4(d). It is seen that the behavior of the 
system in the former (latter) case is periodic (chaotic). Closer to the ignition point, 
the system spends longer time near the lower state and the large amplitude oscillatory 
bursts become less frequent. This and the increase in the size of the attractor across the 
ignition point (occurrence of the large amplitude peaks) can be explained by invoking 
the reinjection property for one-dimensional maps[52, 53]. Here a trajectory leaving the 
vicinity of the lower IP state, like it does when it exhibits a large amplitude oscillation, 
gets reinjected back close to the lower IP state soon thereafter. This principle has been 
discussed in the context of homoclinic tangency and homoclinic orbits by Gaspard and 
Wang[54]. They show how such a behavior can explain the occurrence of mixed-mode 
oscillations. 

The intermittency appears to be the fore-runner of a sequence of more complex 
dynamic behavior which we see in the region 0.415 < P < 0.70. In this region the 
dynamic behavior alternates regularly between periodic and chaotic as the parameter 
P is varied. The periodic states are typical of mixed-mode oscillations consisting of 
a mixture of large amplitude and small amplitude oscillations. In this region we have 
observed a period adding sequence, different period doubling cascades as well as saddle- 
node bifurcations. A review detailing these characteristics can be found in Swinney[49]. 

The characteristics of the stable state of the system in the region 0.415 < P < 0.70 
are summarized in Table 2.3. Each chaotic state is represented by the letter C. Simi- 
larly, each periodic state is represented by the letter P. Each sub-harmonic oscillation 
having a period m times the forcing period is represented with the number m preceding 
the letter P. Each P is followed by one or more subscripts and superscripts. The sub- 
script (superscript) denotes the number of consecutive small (large) amplitude peaks 
in the periodic state. Thus a 11P 9 2 periodic state is characterized by two consecutive 
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Figure 2.4: Variation of the temperature T with dimensionless time ((a) and (b)) 
and the stroboscopic map ((c) and (d)) near the ignition point showing intermittency 
behaviour, parameter set 1, R = 0.573875, (a), (c): P = 0.4105; (b), (d ):' P = 0.411. 
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Table 2.3: Period-adding and alternating periodic-chaotic regions, set i, 11 = 0.573875. 

large amplitude peaks and nine consecutive small amplitude peaks. Fig. 2.5 shows the 
dependence of the largest Lyapunov exponent A on P. Chaotic dynamics are charac- 
terized by a sensitive dependence on initial conditions. This implies that one of the 
Lyapunov exponents is greater than zero. Since the system is dissipative, the other 
Lyapunov exponent is negative. The Lyapunov exponents were calculated using the 
method of Wolf et al. [56]. 

The main characteristic of Table 2.3 is the occurrence of chaotic behavior in between 
regions where the behavior is periodic. This occurrence of alternating periodic-chaotic 
behavior has been observed in several other systems[9, 57-59]. We have observed that 
the chaotic behavior in this region is generated through the period-doubling route. 

The system exhibits a sequence of mP oscillations for successive integer values of 
m (22 to 6 are shown in Table 2.3). The periods of two successive states belonging to 
this basic sequence differ by IP. Each periodic state exhibits 2 large amplitude peaks. 
Transition from one periodic state to another is initiated through a series of period 
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doubling bifurcations to chaos, followed by a recovery to the next periodic state in the 
period-adding sequence[58, 60]. Frequency locking type interaction of two neighbouring 
periodic states yields a new periodic state not belonging to this period-adding sequence. 
Here, between an nP and an (n-l)P state, there is a (2n-l)P state. The (n-l)P and 
(2n-l)P states again interact to give a (3n-2)P state. Repetitive interactions of this 
nature give rise to stable periodic states of different periods. The region of stability of 
these periodic states with large periods is very small[9, 58]. We explain this next in 
the context of Table 2.3. 

For P = 0.56, we see a 15P 6 2 6 state. This is a 15P state with 6 small amplitude 
peaks, followed by 2 large amplitude peaks. This is followed by another sequence of 
6 small amplitude peaks and a large amplitude peak. This particular state arises in 
the region between the 7 P 5 2 and 8P 6 2 states which occur on either side of this value 
of P. We have observed that this 7P state becomes unstable through a limit point 
bifurcation. This unstable 7P state is of the form 7P£. The stable 8P state is of the 
form 8 Pq. It appears that the unstable 7P state interacts with the stable 8P state and 
yields the stable 15P state. Similarly, we observed the existence of a 22P state between 
the 15P and the 7P states. This corresponds to a (3n - 2) periodic state (discussed in 
the preceding paragraph) with n — 8. Similarly, 13P and 19P states exist between the 
7P and the 6P states. We have also found a 25P state between the 13P and the 12P 
states of the period-adding sequence. 

The existence of these periodic states can be explained with the help of the frequency 
locking type of interaction discussed earlier. It appears that for this interaction to 
occur, one of the members of the period-adding sequence must undergo a saddle-node 
bifurcation. We have confirmed this for the 7P state near P = 0.56, in the vicinity of 
which a stable 15P state exists. Similarly, it appears that the 6P and the 12P states 
undergo a saddle-node bifurcation to give stable 13P and 25P states respectively. This 
appears to be consistent with the total number of large and small amplitude peaks 
occurring in the 12P, 13P, 25P and the 6P, 7P, 13P states respectively in Table 2.3. 

Kubickova et al. [21] have obtained bifurcation diagrams for a number of mP states 
as a function of P = R. They have observed saddle-node bifurcations for these periodic 
states. This provides further evidence that such bifurcations can occur for the different 
periodic states in the period-adding sequence, thereby causing an interaction between 
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Figure 2.6: Stroboscopic map of the chaotic attractor, parameter set 1, P = 0.635, 
R = 0.573875. 

two neighbouring periodic states. 

Near P = 0.628, the 7P state undergoes a period doubling bifurcation. This gives 
rise to a stable 14P state. As P is increased, this 14P state again undergoes a period 
doubling bifurcation giving a stable 28P state. The period doubling of different mP 
states of periods m — 8, 9 and 10 have also been observed. It is very likely that there 
are other cascades of period doubling bifurcations of each of the members of the period 
adding sequence. This is another mechanism to generate the chaotic behavior (see 
Table 2.3). The stroboscopic map of a chaotic attractor for P = 0.635 is shown in Fig. 
2.6. This arises from the period doubling of the sub-harmonic 7P state. This figure 
can be viewed as being composed of 7 islets. 

With a view to determining the effect of the parameter R on the dynamic behavior 
of the FBR, we carried out extensive simulations for set 1 for R = 0.2, 0.3 and 0.75. 
The bifurcation diagram depicting the dependence of the IP state for R = 0.20 is also 
S-shaped (Fig. 2.7). The lower branch of the IP state has an ignition point at P — 
0.12051. On the upper branch, the IP state is unstable between 0.1337 < P < 0.20475. 
To the left of this interval, a complex-conjugate pair of eigenvalues leaves the unit circle. 
At the right end point of this interval, a complex-conjugate pair of eigenvalues re-enters 
the unit circle. These points correspond to torus bifurcation points. In this region we 
have observed the existence of quasi-periodic oscillations for 0.1337 < P < 0.1365. 
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Figure 2.7: Bifurcation diagram of the IP state for the FBR, parameter set 1, R = 0.20 
and R = 0.30, ( — ): stable IP state; (•••): unstable IP state; (•): Torus bifurcation 
point. 

The' FBR exhibits periodic, quasi-periodic and chaotic dynamic behavior in the region 
0.1385 < P < 0.20475. Chaotic behavior in this region appears to arise from a periodic 
- quasi-periodic - chaotic transition [49]. Figure 2.8(a) shows the stroboscopic map 
for the quasi-periodic behavior of the FBR with P = 0.13385. As expected, the map 
yields a closed curve. In Figs. 2.8(b) and 2.8(c), we show the stroboscopic maps of 
the dynamics for P = 0.13850 and 0.20 respectively. For P = 0.13850, we see that 
the closed curve has undergone ‘folding’. This is a typical feature of the stroboscopic 
maps in the region 0.1365 < P < 0.1385. As P is increased past this value, the closed 
curve breaks up indicating a transition to chaotic behavior. The mapped points now 
lie on segments of disjoint curves. The onset of chaos is thus signalled by the folding 
and subsequent breaking up of the closed curve. Figure 2.8(c) is representative of the 
typical nature of the stroboscopic map of the behavior obtained for values of P in the 
region 0.1386 < P < 0.20475. Further simulations indicate that the stroboscopic map 
remains similar to Fig. 2.8(c) till the second bifurcation point at P = 0.20475, where 
the IP state regains its stability. 

The bifurcation diagram of the IP state of the FBR for R — 0.75 is shown in 
Fig. 2.9. The qualitative nature of the bifurcation diagram is very similar to that for 
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Figure 2.9: Bifurcation diagram of the IP state for the FBE, parameter set 1, R = 
0.750, (— ): stable IP state; (• • •): unstable IP state; (x): Period doubling bifurcation; 
(•): Torus bifurcation point. 

R = 0.573875. We have a region of multiple IP periodic states for 0.026 < P < 0.5776. 
On the upper branch of the IP state a period doubling cascade begins at P = 0.299. 
This proceeds according to the Feigenbaum scenario. Beyond the accumulation point 
of this cascade, the dynamics are chaotic and the stroboscopic maps are again ’islet’ 
shaped. These islets grow in size and finally merge as P is increased, in much the 
same way as we saw for R = 0.573875 (Fig. 2.3(a)-2.3(c)). Simulations with P beyond 
0.330 always end up on the lower stable IP state. This was tested for a wide variety 
of initial conditions. It appears that the chaotic attractor has disappeared as a result 
of a boundary crisis. As in the case of R — 0.573875, this is due to the occurrence of 
a crisis. This persists till we cross the ignition point ( P = 0.5776). Just beyond the 
ignition point, the behavior appears to be typical of intermittency, i.e., mainly small 
amplitude periodic behavior is interrupted by bursts of large amplitude oscillations. 
It appears that the dynamics of the FBR are essentially periodic in this region. This 
dynamic behavior consists of very large period oscillations. For example, simulations 
have revealed that at P = 0.5777, a 160P state exists. In the region 0.5777 < P < 
1.340, we have observed a period-adding sequence and regions of alternating periodic- 
chaotic behavior. The dynamic behavior is similar to that with R = 0.573875. A 
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complex-conjugate pair of eigenvalues re-enters the unit circle at P = 1.34 and the IP 
state regains its stability. 

The bifurcation diagram of the IP state of the FBR for R = 0.30 is shown in Fig. 
2.7. For this value of R, the dynamic behavior of the FBR in the region on the upper 
branch where the IP state is unstable (0.2164 < P < 0.3295) is similar to that for 
R = 0.20 and R = 0.573875. As P is increased beyond 0.2164, a complex conjugate 
pair of eigenvalues leaves the unit circle. The dynamics of the FBR are quasi-periodic 
for this value of P . With further increase in P, the quasi-periodic dynamics disappear 
and the FBR exhibits chaotic dynamics. With further increase in P beyond 0.220, 
we see a period-adding sequence and regions of alternating periodic-chaotic behavior. 
Cascades of period doubling bifurcations also occur in the period-adding sequence. 

For R — 0.2 and 0.30, the region where the upper IP state branch is unstable occurs 
to the right of the ignition point. Hence in these cases, we do not expect nor observe 
any crisis (Table 2.4). 

Simulations for values of R beyond 0.90 have shown a decrease in the region where 
multiple IP states exist. There is also a decrease in the size of the region on the upper 
branch where the IP state is unstable. This conforms to the fact that the dynamics of 
the FBR approach that of the batch reactor for high R. For set 1, the high value of B 
perhaps necessitates going to such high values of R, in order for the transition to the 
batch reactor to be realized. 

For the parameters in set 2 of Table 2.1, the CSTR has a region of multiple steady- 
states and a region where the upper steady-state branch is unstable and has limit 
cycles [6]. The bifurcation diagram for this set of parameters is shown in Fig. 2.1. The 
dependence of the IP states of the FBR on P has been obtained for three values of R, 
viz., 0.25, 0.50, 0.75. This is shown in Fig. 2.10. 

For R = 0.250, we have a region of multiple IP behavior in the interval 0.0249 < 
P < 0.2715. In the upper branch the stable IP behavior becomes unstable via a torus 
bifurcation at P = 0.227. At P = 0.3103, a complex conjugate pair of eigenvalues 
re-enters the unit circle and the IP state becomes stable again beyond this value of 
P. In the region 0.2715 < P < 0.3103, immediately to the right of the ignition point, 
simulations reveal the existence of quasi-periodic behavior for several values of P. A 
typical stroboscopic map of such behavior for P = 0.28 is shown in Fig. 2.11. The 
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Set 

R 

Behavior observed 

1 

0.20 

Multiplicity, periodic, quasi-periodic and chaotic behavior 


0.30 

Multiplicity, quasi-periodic, chaotic behavior, period-adding 
giving rise to alternating periodic-chaotic behavior 


0.573875 

Multiplicity, period doubling cascade to chaos, crisis, intermittency, 
mixed-mode oscillations, period-adding giving rise to alternating 
periodic-chaotic behavior 


0.75 

Multiplicity, period doubling cascade to chaos, crisis, intermittency, 
mixed-mode oscillations, period-adding giving rise to alternating 
periodic-chaotic behavior, 

2 

0.25 

Multiplicity and quasi-periodic behvaior 


0.50 

Multiplicity and quasi-periodic behavior 


0.75 

Unique stable IP state 

3 

0.25 

Multiplicity and isola 


0.50 

Multiplicity and isola 


0.75 

Unique stable IP state 


* 


Table 2.4: FBR bifurcation behavior. 



Figure 2.10: Bifurcation diagram of the FBR for parameter set 2, (— ): stable IP state; 
(•••): unstable IP state; (•): Torus bifurcation point. 
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Figure 2.11: Stroboscopic map for the FBR showing quasi-periodic behavior, parameter 
set 2, P = 0.28, R = 0.25. 

points trace out a continuous closed curve . In the region 0.227 < P < 0.2715, 
immediately after the IP state becomes unstable on the upper branch, s imula tio ns 
with different initial conditions always end up on the lower stable IP state. This may 
be because the quasi-periodic attractor has been destroyed in a collision with the IP 
periodic state at the saddle-node bifurcation point. 

We notice a similar qualitative behavior for R = 0.50. For R = 0.750, the region of 
multiple IP states and the region where the IP state is unstable disappears (Fig. 2.10). 
The dynamics of the FBR, therefore approach the dynamics of the batch reactor for 
this and higher values of R. Unlike in set 1, we have not observed any period doubling 
bifurcations on the upper IP state branch of the bifurcation diagram for all the values 
of R studied. 

For the par am eter values in set 3 of Table 2.1, the CSTR bifurcation diagram shows 
a region of multiple steady-states and an isola[6]. The bifurcation diagram for this set 
of parameters is shown in Fig. 2.12. For the FBR, we get a similar dependence of 
the basic IP state on P for R = 0.10 and 0.20. The dependence of the IP behavior of 
the FBR for R = 0.1, R = 0.20 and R = 0.50 is shown in Fig. 2.13. For the first 
two values of R, we see a region of multiple IP states and an isola. As we increase 
R from 0.10 to 0.20, the regions of existence of multiple IP states and of the isola 
decreases. When R is increased further to 0.50, the isola and the multiple IP states 
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Figure 2.12: Bifurcation diagram of the CSTE for parameter set 3, (— ): stable steady- 
state; (•••): unstable steady-state; (x): Hopf bifurcation point. 



Figure 2.13: Bifurcation diagram of the FBR for parameter set 3. ( — ): stable IP 
state; (• • •): unstable IP state. 
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disappear. A unique stable IP state exists for all values of P. Thus, again we approach 
the dynamics of the batch reactor at high values of R. Table 2.4 contains a summary 
of the bifurcation features we have observed for the FBR. 

2.5 Conclusions 

In this work we have presented bifurcation features of the FBR for selected values of 
the parameters B, fio and /cq. Our objective was to study how the local bifurcations in 
the CSTR manifest themselves in the FBR. From the numerical evidence presented in 
the previous section, we may conclude that the steady-state bifurcation (saddle-node 
bifurcation) in the CSTR manifests itself as a saddle-node bifurcation in the FBR. The 
Hopf bifurcation in the CSTR can manifest itself as either a period doubling bifurcation 
or a torus bifurcation. It appears that complex, i.e., sub-harmonic, quasi-periodic and 
chaotic, oscillations occur in the FBR only in cases where the equivalent CSTR has a 
Hopf bifurcation. It also appears that the nature of the steady-states of the CSTR also 
influences the dynamics of the FBR. The node and saddle type steady-states manifest 
themselves in a similar way in the FBR. However the complex dynamics of the FBR 
appear to arise only when the equivalent CSTR has a focus type of steady-state. 

A positive feature of this system is that it exhibits complex dynamic behavior over 
a sufficiently wide range of parameter values. In contrast, Jorgensen and Aris[9] report 
complex bifurcation sequences occurring over very narrow ranges of the bifurcation 
parameter. 

The quasi-periodic, chaotic dynamics and regions of period-adding sequence re- 
ported above were obtained by simulating with a wide variety of initial conditions. 
However, the existence of other attractors in these regions cannot be ruled out[21]. 
The basin of attraction of these other attractors may be very small. A thorough in- 
vestigation of the effect of initial conditions will be necessary to verify the existence of 
other attractors. 

For this reactor system, we have shown the existence of a wide spectrum of complex 
dynamic behavior, i.e., the period doubling route to chaos, period-adding sequences 
resulting in alternating periodic-chaotic behavior and period doubling bifurcations of 
various members of the period-adding sequence, quasi-periodic oscillations, crisis and 
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intermittency behavior (Table 2.4). Such dynamic behavior is a universal character- 
istic of nonlinear systems of varying degrees of complexity, from the one-dimensional 
logistic map to systems of coupled ordinary differential equations. Nonlinear electrical 
circuits and chemical reactions have been known to exhibit some of these features. 
This behavior has been analysed in the literature, mainly in the context of discrete 
dynamical systems. The understanding of these features in the context of continuous 
time systems remains incomplete. 

We have compared a CSTR operating at a particular residence time r with the 
FBR with the same ratio of Pj R. There is an infinity of P and R values of the FBR 
which yield a given r. Hence the operation of a CSTR can be compared with an infinity 
of FBR’s. As explained earlier, the FBR approaches the CSTR in the limit P — > 0, 
R —»• 0 and for finite r. 

The periodically forced CSTR has been studied extensively in the literature. The 
unforced CSTR itself may have limit cycles and the interaction of the natural and 
forcing frequencies gives rise to complex dynamics. We have studied the forced op- 
eration of a batch reactor. The type of forcing used in the present study is discrete. 
This causes a discontinuous variation in the concentration and temperature. We have 
observed that, eventhough there is a basic difference in the way the system is being 
forced, the dynamic behavior is equally rich. It must be noted that that the system we 
are forcing, the batch reactor does not show any oscillations, whereas the CSTR our 
system tends to, can exhibit autonomous limit cycles. 

We can also consider the FBR to be the most general form of the well-stirred reactor 
model. This is because both the batch reactor and the CSTR are the two limiting cases 
of FBR operation. From this point of view, it is perhaps natural that the FBR shows 
such complex behavior. 



Chapter 3 


Steady-state Behavior of a 
Reactor-Separator System 

Nomenclature 

B Dimensionless heat of reaction. 

Dao Damkohler number, reaction rate constant. 

Fq Molar flow rate of fresh feed to reactor. 

F Molar flow rate of feed to flash unit. 

k Equilibrium constant. 

L Molar flow rate of the liquid stream from flash unit. 

M r Molar holdup of the reactor. 

P Pressure in the separator. 

P sat Saturated vapor pressure. 

R Molar flow rate of the recycle stream. 

T Dimensionless temperature. 

Tp Dimensionless temperature of the separator. 

V Molar flow rate of the vapor stream from the separator. 

x e Composition of the liquid (recycle) stream, mole fraction. 

x R Composition of the material in the mixer, mole fraction. 

y e Composition of the vapor (product) stream, mole fraction. 

z Composition of the material in the reactor, mole fraction. 
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Greek letters 

a 

Recycle rate per unit molar flow rate of fresh feed. 

A) 

Heat removal parameter. 

Sx e 

Constant in section 3.3. 

P 

Molar density of mixture. 

T 

Mr/Fi o- 

T f 

m r /f. 

n 

m r /l. 

Subscripts 

i 

Component i. 

0 

Conditions at reactor inlet. 

1 

Reactant. 

2 

Product. 


3.1 Introduction 

In the previous chapter we investigated the dynamic behavior of a discretely forced 
reactor system. Another class of systems which have only attracted attention in recent 
years are the reactor-separator systems coupled by recycle. In a continuing series of 
papers, Luyben and co-workers have investigated the steady-state features of several 
simple reactor-separator systems. However, as we pointed out in the first chapter, 
they have not made any attempts to carry out a systematic study and generalize 
the qualitative features of such systems. Recently Morud and Skogestad[25, 26] have 
presented reviews of the effect of recycle on the dynamics and control of plants. With 
the help of several examples, they catalog the different features of linear recycle systems. 
They also highlight the need for a systematic classification and analysis of effects which 
may introduce complex behavior in a plant. This is essential in order to enhance our 
understanding of the different phenomena and also to provide some indicators of when 
the phenomena occur. 

In a series of papers, Luyben and co-workers [28-30] have studied the dynamics 
and control of binary and ternary recycle systems. They consider a process consisting 
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of an isothermal CSTR coupled through recycle to a distillation column. The effect 
of design on the steady-state economics and the controllability of the process were 
discussed. Starting with a given set of design parameters, their approach consists of 
developing a design procedure for the operation of the distillation col umn . The effect 
of changing the design parameters on the system variables was then studied. Based 
on their results, they have conjectured about the advantages of working with certain 
sets of design parameters. However, the complexity of the system in ter ms of the large 
number of variables, makes it impossible to obtain a generalization of their results. 
In a later paper[31], they analysed the effect of different control configurations on the 
snowball effect. The snowball effect manifests itself as a very large change in the 
recycle flow rate for a small change in the independent variable. However, they did not 
provide any explanation about the origin of this phenomenon and ways for preventing 
its occurrence. Recently, they have extended their investigations to cover different 
cases of ternary systems [32-34]. In these studies they consider several model reaction 
schemes, multiple separator units and multiple recycle streams. 

In this chapter we study the steady-state characteristics of a coupled reactor- 
separator system. The reactor is assumed to be a non-adiabatic CSTR. It sustains 
the exothermic, irreversible first order reaction A-* B. The mixture leaving the reac- 
tor is sent to a separator which is modelled as a single stage flash unit. The reactor 
effluent stream is separated into a reactant rich liquid fraction and a product rich vapor 
fraction. The coupling between the two units is achieved through the recycle of the 
reactant rich liquid stream from the separator to the reactor (Fig. 3.1). 

We study the behavior of this model nonlinear reactor-separator system. We de- 
termine the qualitative features of the steady-state behavior of the system for different 
possible combinations of the specified variables. This enables us to identify different 
combinations of the specified control variables which permit operation of the coupled 
system with separation, at a steady-state. A comprehensive analysis is carried out for 
the following sets of operation of the separator 

1. Constant temperature and pressure, (T F , P) flash, 

2. Constant temperature and constant liquid recycle rate, (Tp, L) flash. 

In the first case, operation of the flash is isothermal and isobaric. The vapor-liquid 
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equilibrium (VLE) considerations for the binary system then uniquely determine the 
compositions of the liquid and vapor streams. For this choice of the specified variables, 
separation into a product rich vapor fraction and a reactant rich liquid fraction, is 
possible only when the reactor effluent composition lies in between the equilibrium 
compositions of the liquid and vapor streams. As we shall show later, there exist 
choices of the specified variables for which this constraint is not satisfied. This results 
in the existence of regions in parameter space for which steady-state operation with 
separation is not possible. 

In the second case, the flow rate of the liquid stream from the flash is fixed. Here, 
the pressure in the separator is not externally imposed. Consequently, the compositions 
of the overhead and bottom streams leaving the flash are not fixed. Since the recycle 
stream flow rate L is fixed, the system is forced to regulate itself such that the recycle 
rate is always at its specified value. This ensures that steady-state operation with 
separation is possible for all parameter values. 

For each of the above cases, we study the different possible steady-state multiplicity 
behavior of the coupled system using singularity theory. Application of singularity 
theory enables us to classify a suitable parameter space into regions having different 
multiplicity features. The results are presented in the form of bifurcation diagrams. 
We also make a qualitative analysis of the dynamic behavior which can be expected 
once a choice of the set of variables to be specified has been made. 

The model is described in the next section. In section 3.3 we study the steady-state 
multiplicity behavior for different sets of specifications of this coupled system, as a 
single control parameter is varied. 

3.2 Description of the Model 

The system under consideration is shown in Fig. 3.1. Fresh feed containing pure A 
is fed at a rate Fo to the reactor. The effluent stream from the reactor passes through 
a heat exchanger which changes the temperature of that stream to the specified flash 
temperature Tp. This stream is then flashed into two fractions. The liquid stream 
is rich in the reactant, while the vapor stream is rich in the product. The vapor 
stream is drawn off and the liquid stream is recycled to the reactor through a mixer 
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Figure 3.1: Schematic diagram of the reactor-separator system with recycle. 

(where it mixes with the liquid already present) and a heat exchanger, which brings 
the temperature of the recycle stream to that of the fresh feed. The presence of the 
heat exchangers ensures that there is no thermal coupling between the two units. For 
simplicity, we neglect the dynamics of the heat exchangers. The coupling between the 
two units is therefore, only through the composition and the flow rate of the recycle 
stream. 

Throughout this study, we assume the molar holdups of the reactor and the mixer 
to be constant. The material balance and the equilibrium relations modelling the 
separator usually involve the species composition expressed in terms of mole fractions 
The mass balance in the reactor is normally expressed in terms of the concentration 
of the different species a. To represent both the units on a common basis, we use mole 
fractions to model the system. This exploits the relationship 


We assume that the molar density in the liquid phase p is constant (i.e., independent of 
composition). The overall material balance and the dimensionless material (mole) and 
energy balance equations governing the evolution of the composition of the reactant a, 
and dimensionless temperature T in the reactor are given by the following equations 
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respectively 


F 

= Fo + R 

(3.2) 

dz i 
dt 

R T 

= 1 -Zi - — {zi -x R ) -Da Q re T Zi 

(3-3) 

dT 

dt 

= — T(1 + + M + BDa Q re T zi 

-CO 

(3-4) 


Here B, fa and Da 0 represent the dimensionless heat of reaction, the dimensionless 
heat transfer coefficient and the Damkohler number (representing the reaction rate 
constant) respectively. We have used r = Mr/ Fq as the characteristic time to define 
the dimensionless time. We have made the exponential approximation in the derivation 
of the above equations. The overall and reactant mole balances for the equilibrium 
separation process are given by 

F = V + L (3-5) 

Fz x = Vy e i + Lx e i (3-6) 

Subscript 1 (2) refers to the reactant (product). F, V and L represent molar flow rates 
of feed, vapor and liquid streams to and from the flash. The reactant rich liquid stream 
passes through a constant holdup mixer before being recycled to the reactor. Therefore, 
the recycle rate to the reactor is equal to the rate of the liquid stream entering the 
mixer (i.e., R = L). However, this rate is not constant but depends on x eX , y ei and z x . 
The evolution of the reactant composition in the mixer is governed by the following 
dimensionless equation 

IT = “ xm) (3 ' 7) 

The separation process is governed by the vapor-liquid equilibrium relations relating 
the compositions (x e i, y e i, i = 1,2) in both the phases. We assume the vapor and 
liquid phases to be ideal, so that Raoult’s law is valid. The sum of the mole fractions 
in the liquid and the vapor phase must equal unity. Therefore, we have the following 
additional sets of equations 


Vel — k\X & i 
Ve 2 = ^2^e2 


(3.8) 

(3.9) 
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2/el + Ve2 = 1 

x el + x e2 = 1 


(3.10) 

(3.11) 


ki and k% are the equilibrium constants defined as ki = Pf at /P, where P? at is the 
saturation vapor pressure of component i. The saturation vapor pressures are assumed 
to be given by the Antoine relations. Then using (3.2), (3.5) and (3.6), the recycle rate 
per unit fresh feed entering the reactor, a, can be defined as 

R_ _ Vel ~ Zl 

F 0 ~ 


L 

a = — 
F 0 


(3.12) 


Zi X e \ 

The set of equations (3.2)-(3.11) are a set of 10 equations in the 14 unknowns: Fq, F, 
L, P, T, Tp , V, Mr, x e i, x e2, x r, Vei, Ue2 and Z\. Hence four additional variables need 
to be specified so that the above equations can be solved for the remaining variables. 

In the following analysis, all the compositions are expressed as mole fraction of the 
reactant A, and hence the subscript 1 is suppressed. 


3.3 Steady-State Analysis 

In the first part of this section we consider the isothermal and isobaric {Tp, P) operation 
of the separator. This choice fixes two variables. We consider the following sets of 
specifications for specifying the remaining two variables. 

1 . {M r , F 0 ) fixed: The reactor holdup Mr and the fresh feed flow rate F 0 axe spec- 
ified, 

2. {Mr, F ) fixed: The reactor holdup Mr and the reactor effluent stream flow rate 
F are specified, 

3. {Mr, L) fixed: The reactor holdup Mr and the recycle stream flow rate L are 
specified. 

4. {F 0 ,F) fixed: The fresh feed flow rate F 0 and the reactor effluent stream flow 
rate F are specified. 

The overhead (vapor) and bottom (liquid) stream compositions for the case of isother- 
mal and isobaric separation can be determined using (3.8)-(3.11). The model equations 
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(3.2)-(3.11) can then be reduced to a system of three coupled ordinary differential equa- 
tions in z,T and xr. The resulting model is solved for the steady state by setting the 
time derivatives to zero. 

In subsection 3.3.5 we consider the flash temperature Tp and the liquid stream 
flow rate from the separator L to be fixed. We also consider the molar holdup of 
the reactor Mr and the fresh feed flow rate To to the reactor to be fixed. With this 
set of specifications, the pressure of the separator remains an unknown quantity. The 
equilibrium constants ki and k 2 are not known and the equilibrium compositions x ei 
and y e i have to be calculated by solving the set of equations (3.2)-(3.11) simultaneously. 

3.3.1 Isothermal and Isobaric Separation with (Mr, Fo) Fixed 

We begin by considering the case when the holdup in the reactor Mr and the fresh feed 
flow rate Fq are the specified variables. The recycle stream flow rate and the reactor 
effluent stream flow rates are now unknown and are dependent on the composition 
in the reactor. At steady-state, the resulting system of coupled nonlinear algebraic 
equations (3.2)-(3.11) can be reduced to the following single equation in the unknown 
T 


H(T,p , r) = Daor(x e — y e )Te T + x e Daorl3orTe T — x e BDa Q r(l — y e )e T 

—Por(l — y e )T + B(l — y e ) 2 = 0 (3.13) 


In the unknown z, the corresponding equation is 


G(z,p,r ) = In 


' 1-2/e ' 

.DoqTZ. 


Xe -Ve + A) r{x e - z)} - B( 1 - y e )(x e - z) = 0 


(3.14) 


Here p represents the parameter set ( B,j3 0 ,Da 0 ,x e ,y e ). The parameter r is defined 
as the ratio of the molar holdup in the reactor Mr and the molar flow rate of fresh feed 
To to the reactor, r is considered to be the bifurcation parameter and is hence given a 
special identity in (3.13) and (3.14). The residence time in the reactor is dependent on 
To and the recycle rate R, which in turn is dependent on the reactor effluent stream 
composition. It must be emphasised that r does not represent the residence time in 

the reactor, since it depends only on Fq and not on the reactor effluent stream flow 
rate. 
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The nonlinear equations (3.13) and (3.14) in the unknowns T and z respectively can 
have multiple solutions for a fixed set of parameters. We investigate the multiplicity 
features of (3.13) and (3.14) using singularity theory[35]. Singularity theory enables us 
to obtain a complete picture of the different possible steady-state multiplicity features 
of a system. A brief summary of singularity theory is presented in the appendix. We 
now present the singularity theory analysis of (3.13) and (3.14) in detail. 

The set of equations (the subscripts T, z and r indicate differentiation w.r.t. the 
respective variables) 

H = H t = Htt = 0 (3.15) 

has the solution 


N 


4- 


4x e B{l - y e ) 


r = 


B{ 1 - y e ) 


(3.16) 


(%e Ve) Po(T 2 ) 

This solution is infeasible, since T < 2 and hence r < 0. The corresponding solution 
for z satisfying the set of equations 


G = G Z = G, 


0 


(3.17) 


is given by: z = (1 — y e )/DaoTe T , with T and r given by (3.16). This implies that 
(3.13) or (3.14) can only have a maximum of two solutions for any given set of the 
operating parameters. Equations (3.13) and (3.14) behave like quadratic equations, 
and therefore, there exists the possibility that for certain operating conditions, the 
system will not have any steady states. This is discussed in detail later. Such behavior 
is in contrast to the CSTR where we can have a maximum of three solutions for any 
given set of parameters[6, 23]. 

Further, the set of equations 


G = G z = G r = det[d 2 G\ = 0 

G zz G. 


where det[d?G] is the determinant of the matrix 
has the solution given by 


ZZ ZT 

G rz Gtt 


x. 


y e = z = 


Po 


BDa oe 1 


, T = 1, r = 


B(l- Ve) 


(3.18) 


(3.19) 



After specifying B, Pa and Da 0 , we solved (3.18) for the unknowns z, r, x e and y e . 
The solution to (3.18) could not be obtained analytically. The numerically computed 
values of x e and y e satisfying this set of equations appeared to be equal to each other 
for several values of B, Pa and Daa- It was possible to solve and obtain the singular 
point (3.19) after substituting x e = y e in (3.18). This singular point is a point in 
parameter space p where the function (3.14) is degenerate, i.e., it satisfies the equality 
conditions (3.18). The system has a highest order singularity when x e and y e take the 
above values for fixed values of the parameters [B, Pa, Daa). 

Irrespective of conditions in the reactor, the compositions of the vapor and liquid 
streams exiting the flash are fixed, since we have assumed isothermal and isobaric 
separation. A steady-state is defined to be feasible if it permits steady-state operation 
of the coupled system with separation. Hence a steady-state of the coupled system is 
feasible only ii y e < z < x e . Therefore, the singular point (3.19) lies in the infeasible 
region of operation, i.e., it cannot be achieved in practice. 

When x e = y e = z, the separator acts as a stream splitter. The solution of (3.14) 
given by (3.19) does not model the physical process under consideration (reaction 
and separation). Under these conditions, Equations (3.5), (3.6) and (3.8)-(3.11) are 
no longer valid, as no separation is effected. The reactor system can still be modelled 
using (3.3) and (3.4). However, we need to fix (arbitrarily) R/Fa ■ The model equations 
(3.3) and (3.4) would now represent the dimensionless material and energy balance 
equations for the CSTR without recycle (R/Fq = 0) and with recycle (for a positive 
non-zero R/Fa). The implications of the presence of infeasible regions on the steady- 
state behavior of the coupled system are discussed next. 

As the steady-state value of z approaches y e , more and more of the feed to the 
separator goes into the vapor phase and the recycle rate approaches zero. Consequently, 
in the limit of z = y e , the coupled system reduces to a system connected in series. In 
this case, the flash unit does not affect the behavior of the CSTR. Further, as the 
steady-state value of z approaches x e , i.e., for low conversions in the reactor, a larger 
proportion of the feed entering the flash goes into the liquid phase and the recycle rate 
(Equation (3.12)) tends to infinity. This results in a paradoxical situation in the limit 
°f z — x e, where we have total recycle of the reactor effluent stream. This contradicts 
the assumption of steady-state operation under constant molar holdup of the reactor 
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and the mixer, with fresh feed also entering the reactor at a constant rate. To overcome 
this situation, we impose the following bounds on the conversion in the reactor, for a 
feasible steady-state operation of the coupled system 

y e < z < x e - Sx e 

Here 6x e is a small positive constant incorporated to exclude the possibility of total 
recycle in the process and its governing system of equations. This bound is imposed 
on the conversion to enable a feasible steady-state operation of the coupled system. 
The inclusion of a non-zero Sx e is equivalent to specifying an upper bound on the 
recycle stream flow rate. An additional feature of singularity theory is that it can 
systematically account for the changes in the multiplicity features resulting from the 
presence of such feasibility boundaries on the state variables. 

From (3.12), it can be seen that in the limit z x e , a small change in the fresh feed 
stream flow rate is associated with a large change in the recycle stream flow rate. This 
is the snowball effect observed by Luyben[30]. In general, the snowballing increases in 
magnitude only when perturbations acting on the system cause 2 to approach x e . 

To determine the normal form of the singularity, we evaluate the eigenvector v 
corresponding to the zero eigenvalue of the Hessian matrix [d?G]. The only non-zero 
element of the Hessian matrix [d 2 G] is G zz = 2 (5 0 r/z. The eigenvector v is therefore, 
determined to be [0 1] T . It is also necessary to determine the directional derivatives of 
G along u[35]. The directional derivative G v is defined as 


G v = G vl = [G z G r ] [0 l] r (3.20) 

At the singular point, G v — 0. Similarly, G vv , G vvv , G vvvv and G vvz are given by 

G vv = G v 2 = [G v iz G„i t ] [0 l] r (3-21) 

G m =G v 3 = [G v2z G v 2r][0 1] r (3.22) 

G vvvv = G v 4 = [G v 3 Z G v zr] [0 l] r (3.23) 

G VV z = [G V 2z G V 2z] [0 l] r (3-24) 


In the above expressions, the subscripts 2 and r indicate differentiation with respect 
to these variables. At the singular point given by (3.19) 

G vv = G V w = Gwvv — 0 5 G V vz — A)/ r an d ( G VVVV G ZZ ) — G vvz ^ 0 (3.25) 
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The normal form of the singularity is x 2 — A 4 = 0[35], where A is the bifurcation 
parameter. This singularity is of co-dimension three and its universal unfolding is given 
by 

x 2 - A 4 + a + f}\ + 7 A 2 = 0 (3.26) 

where a, j3 and 7 are auxiliary parameters. In a neighborhood of the singular point 
given by (3.19), the multiplicity features of the coupled system are similar to the 
multiplicity features of this polynomial equation near the origin in the (a, /3, gamma ) 
parameter space. As mentioned earlier, the singular point (3.19) lies in the infeasible 
region of operation. However, we can still obtain the multiplicity features in a neighbor- 
hood of this singular point. As explained in the appendix, singularity theory enables 
us to identify and locate the critical surfaces surfaces in parameter space across which 
the multiplicity features change[23, 24]. For the set of the specified variables under 
consideration, the hysteresis variety does not exist in the feasible region of parameter 
space (since on this r < 0). Hence the multiplicity features can change only when 
the isola variety is crossed. The nature of the bifurcation across the isola variety is 
determined by the sign of the determinant of the Hessian matrix [d?G\. For positive 
(negative) values of the determinant, we have a simple (isola) bifurcation. Computing 
the isola variety enables us to demarcate a suitable subset of the parameter space p 
into regions having different multiplicity features. As mentioned earlier, the singular 
point of the coupled system is of co-dimension three. For a singularity of co-dimension 
three, all the multiplicity features can only be obtained by representing the varieties in 
a three-dimensional parameter space. However, for ease of representation, the surfaces 
of the different boundaries are projected onto the x e — y e plane. This is equivalent to 
representing the behavior in the (2>, P) plane of flash operation. 

As mentioned earlier, for an isothermal and isobaric flash, a steady-state of the 
coupled system given by (3.14) is feasible only when y e < z < x e — Sx e . For the sake 
of convenience, we take Sx e = 0. These constraints appear as feasibility boundaries 
in the bifurcation diagrams. The steady-states of the coupled system are not feasible 
when z lies outside the feasibility region. New types of multiplicity features arise 
depending on the manner in which the boundaries are crossed. Singularity theory can 
also be employed to account for these new kinds of multiplicity features in a systematic 
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way. Special parameter sets have been defined to account for the presence of the 
feasibility boundaries[24, 35]. These enable us to demarcate the locus of the feasibility 
boundaries in parameter space. These are summarized in the appendix. Next we study 
the existence of these sets for the system. 

• Boundary limit set (BLS): The parameter set corresponding to which a turn- 
ing point of the steady-state branch in the bifurcation diagram lies on the feasi- 
bility boundary forms the BLS. 

For the turning point on the x e boundary, the value of x e is determined to be 

= 2/e/(l -S(l Ve)) (3.27) 

For all y e in the feasible region (0, 1), x e predicted by (3.27) is always outside the 
feasible region. Consequently, the BLS exists only for z = y e . 

• Boundary tangent set (BTS): For parameter values belonging to this set, 
the steady-state branch in the bifurcation diagram is tangential to the feasibility 
boundary. 

The BTS for z = x e corresponds to the diagonal x e = y e . This lies in the 
infeasible region in the parameter space. Whether or not the locus of the BTS 
corresponding to the crossing of the z = y e boundary exists, depends on the 
values of the parameters B, Pq and Da^. 

• Double cross set (DCS): This is given by the values of x e - y e for which two 
branches of the bifurcation diagram cross the feasibility boundaries simultane- 
ously. 

The three boundary sets were computed numerically. The isola variety and these 
three boundary sets divide the x e - y e plane into different regions. For the case of 
(M r , F 0 ) fixed and isothermal and isobaric separation, these are the only surfaces in 
parameter space across which the multiplicity features change. These are depicted in 
Fig. 3.2 and 3.3 for two sets of the parameters B, Po and Da Q . In Fig. 3.2 and 3.3, 
an isola bifurcation occurs along the curve AB, while a simple bifurcation occurs along 
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Figure 3.2: Global classification of the steady-state bifurcation diagrams for (Tf, P) 
and ( M r ,F 0 ) flash (Equation (3.14)) in the x e - y e plane for B = 15, /3 0 = 4.0, 
Da o = 0.20. BLS: Boundary limit set; BTS: Boundary tangent set; DCS: Double cross 
set; IV: Isola variety. 
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BC. The isola variety and the three boundary sets divide the x e - y e plane of Fig. 3.2 
[Fig. 3.3] into twelve (eleven) regions. 

The bifurcation diagrams in the different regions of Figs. 3.2 and 3.3 are shown 
in Figs. 3.4 and 3.5. For a fixed set of the operating parameters, either zero or 
two solutions exist for a given value of r in the bifurcation diagrams. However, the 
feasibility limits on z can give rise to regions where a unique feasible solution exists. 
The changes in the bifurcation diagrams due to crossing of the various boundaries in 
Fig. 3.2 is explained next for a few cases. As we move from region (e) to region (f), we 
cross the BLS (z = y e ). Thus, in the two bifurcation diagrams corresponding to these 
regions [Fig. 3.4(e) and 3.4(f)], the turning point lies above and below the feasibility 
boundary for z = y e . Similarly, as we move from region (b) to region (c), we cross the 
DCS. In both regions, the bifurcation diagrams are such that the number of feasible 
solutions changes as 1 — 2 — 1 as r is varied [Fig. 3.4(b) and 3.4(c)]. The solution 
branch crosses the z = y e boundary first (at a lower value of r) in region (b), while it 
crosses the z = x e boundary first in region (c). 

The bifurcation diagrams on either side of the BTS of Fig. 3.2 are shown in Fig. 
3.4(a) and 3.4(b). In Fig. 3.4(a) the lower steady-state branch lies in the infeasible 
region (z < y e ). A small part of the lower steady-state branch lies in the feasible region 
in Fig. 3.4(b). This is not clearly seen in the figure. When we cross curve AB of the 
isola variety in Fig. 3.2, an isola is born or disappears. This can be seen from the 
bifurcation diagrams in Fig. 3.4(d) and 3.4(e) respectively. In Fig. 3.4(e), there is 
an additional isolated branch in the bifurcation diagram. The curve BC of the isola 
variety represents the parameter set for which a simple bifurcation occurs at a certain 
value of r. The bifurcation diagrams on either side of this curve are shown in Fig. 
3.4(c) and 3.4(d). 

The BLS and BTS corresponding to z = x e - 6x e depicted in Figs. 3.2 and 3.3 
have been obtained by assuming 6x e = 0, i.e., for z = x e . For small nonzero 5x e , these 
sets are shifted slightly from their present position (the diagonal x = y), thus, giving 
rise to new sets of multiplicity features. These new multiplicity features occur in very 
small regions in parameter space and hence are not considered in this study. 

In the case of the CSTR, we are assured of the existence of atleast one steady- 
state for all combinations of parameters. In complete contrast, it can be seen from the 
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Figure 3.4: Schematic bifurcation diagrams of the steady-state for (2>, P) and (M R ,F 0 ) 
flash (Equation (3.14)) describing the dependence of composition in the reactor on r. 
Letters indicate regions in Fig. 3.2. The horizontal lines represent the feasibility 
boundaries z = x e and z — y e respectively (x e > y e ). 
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Figure 3.5: Schematic bifurcation diagrams of the steady-state for (7>, P) and ( Mr,F 0 ) 
flash (Equation (3.14)) describing the dependence of composition in the reactor on r. 
Letters indicate regions in Fig. 3.3. The horizontal lines represent the feasibility 
boundaries z = x e and z — y e respectively ( x e > y e ). 
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bifurcation diagrams that there exist regions of r where no solution exists. This is a 
result of the particular choice of the specified variables. Specifying M R and F 0 and 
having an isothermal-isobaric flash operation constrains the performance of the coupled 
system. This is because the vapor flow rate and vapor composition are determined 
apriori irrespective of the residence time in the reactor. The overall mole balance 
determines V to be equal to F 0 . The VLE considerations determine y e for the binary 
system. As a result, the performance of the coupled system is not allowed to be 
influenced by the interactions between the two individual units. The system behavior 
here is determined solely by the parameters (Tp, P ) and (Mr,F 0 ) and not by the other 
conditions in the reactor. Our results indicate that under some operating conditions, 
such a constrained system may not have any steady-states. 

Luyben[30] has studied the effect of different control structures on the steady-state 
charcteristics of several recycle systems. They report the possibility of occurrence of 
the snowball effect with some control structures. The analysis presented earlier in this 
section, indicates that snowballing increases in magnitude a s z x e . This occurs only 
because x e is fixed, since we assume an isothermal and isobaric flash operation. In 
addition, we have assumed a fresh feed stream entering the reactor at a constant rate. 

The region in the ( x e , y e ) plane that is important in practical situations is the one 
corresponding to x e -» 1 and y e — 0. This corresponds to operation in region (c) and 
(d) of Fig. 3.2 and region (c) of Fig. 3.3. The corresponding bifurcation diagrams are 
shown in Fig. 3.4 and Fig. 3.5 respectively. Operating the system around the lower 
steady-state branch of the bifurcation diagrams should make it possible to avoid the 
possibility of very large internal flow rates, associated with the presence of the snowball 
effect. 

The steady-state equation for the temperature in the reactor (3.13) can be similarly 
analysed. It can be shown that the set of equations 

H = H t = H T = det[d 2 H] = 0 and H vvv ± 0 (3.28) 

also has the solution given by (3.19). The equation (3.13) is contact equivalent to 
x 2 — A 3 = 0 at the singular point. The differences in the normal forms for temperature 
and composition in reactor problems is well known [23]. We have however, not employed 
this equation to determine the multiplicity features as the feasibility constraints do not 
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affect this variable directly. 


3.3.2 Isothermal and Isobaric Separation with (Mr, F) Fixed 


We now consider the case when the reactor effluent flow rate F is the specified variable, 
in addition to the molar holdup in the reactor Mr. The fresh feed flow rate F Q to the 
reactor is now an unknown quantity and must be computed from the model equations. 
Since the flash is isothermal and isobaric, the compositions of the streams leaving the 
flash are fixed. The reactor composition z and temperature T, the recycle stream flow 
rate R and the fresh feed flow rate Fq all change such that the reactor effluent stream 
flow rate F is at its specified value. The dimensionless time in the model governing 
the dynamic behavior of the system (obtained from Equations (3.2)-(3.11)) is now 
defined with respect to the actual residence time in the reactor. At steady-state, these 
equations can be rearranged to yield the following single equation in z 


Gi (z,p,r f ) = B{1 - y e ){z - x e ) - (1 + ftr f)(y e — x e ) 


(1 - y e )(z - Xe) 
Da 0 r f z(y e - x e ) 


(3.29) 


r f in the above equation represents the actual residence in the reactor. The parameter 
set p represents the parameters ( B , ft, Da 0 ,x e , y e ). The solution to the set of equations 


G\ = G u = G i* z = G\ Tf — 0 


is given by 


(3.30) 


Z = X e /2, Tf ~ 


jzVe 

Da 0 (x e - y e )e 2 ’ 


R _ 8 ( a; e ~ Ve) n _ Da 0 (x e - y c )e 2 
x e (l-y e Y 1 -y e 


It can be shown analytically that there is no solution to the winged cusp singularity, 
i.e., there is no solution satisfying the set of equations 


Gi — G\ z — G\ zz — G\ T j — G\ ZTj = 0 (3.31) 

Therefore, for this set of specified variables, the system has a pitchfork singularity at 
its most singular point. 

A typical projection of the numerically computed hysteresis and isola varieties and 
the boundary sets (BLS and BTS for z = y e ) in the (B, ft) plane is shown in Fig. 3.6. 
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Figure 3.6: Global classification of the steady-state bifurcation diagrams for (Tp, P ) 
and (Mr,F) flash (Equation (3.29)) in the B — /3 0 plane for x e = 0.95, y e = 0.1, 
Da o = 0.10. BLS: Boundary limit set; BTS: Boundary tangent set; HV: Hysteresis 
variety; IV: Isola variety; LCS: Limit and cross set. 
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Figure 3.7: Schematic bifurcation diagrams of the steady-state (Tp, P) and (Mr, F) 
flash (Equation (3.29)) describing the dependence of composition in the reactor on 
Tf. Letters indicate regions in Fig. 3.6. The horizontal lines represent the feasibility 
boundaries z = x e and z = y e respectively (x e > y e ) . 
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Equation (3.29) behaves as a cubic equation, and therefore we will have atleast one 
steady-state for all values of the operating parameters, although this can be infeasi- 
ble. For this set of the specified variables, the absence of feasible steady states is a 
consequence of the presence of the feasibility bounds on the conversion in the reactor. 

The coupled system cannot have a feasibility boundary at z = x e . This is because 
fixing the reactor effluent stream flow rate implies that the reactor composition z cannot 
equal x e , since this implies complete recycle of the recycle effluent stream. Complete 
recycle of the product from the separator and the assumption of constant holdup in 
the reactor, force F 0 to be equal to zero. As a result the product stream flow rate also 
equals zero, resulting in the process operating in the batch mode. This in turn leads 
to an increase in the conversion, i.e., a decrease in 2 . 

The coupled system can still behave as a system without recycle, eventhough the 
reactor effluent stream rate has been specified. This corresponds to the case where 
we have complete vaporisation of the feed to the separator. The reactor conversion 
z in this case would then be less than y e . Therefore, in this case only the feasibility 
boundary corresponding to z = y e , exists. 

The bifurcation diagrams in the different regions are shown in Fig. 3.7. From Fig. 
3.7 it can be seen that for some ranges of r/, there are no feasible steady states for 
the coupled system. The multiplicity features of the coupled system are essentially 
similar to that of a CSTR[23]. The isola and hysteresis varieties (Fig. 3.6) reduce 
to those of the CSTR in the limit of x e = 1 and y e = 0. The nature of the changes 
in the bifurcation diagrams when the boundary sets are crossed are explained in the 
appendix. 

Practical considerations for an efficient operation of the coupled system require 
that the feasibility boundaries should be avoided. Prom Fig. 3.6 and the bifurcation 
diagrams shown in Fig. 3.7, it can be seen that this is possible provided we operate 
the system in regions (a), (b) and (j) of Fig. 3.6. 

3.3.3 Isothermal and Isobaric Separation with (Mr, L) Fixed 

We now consider operation of the coupled system when the liquid rate L leaving the 
flash is the specified variable, in addition to the molar holdup in the reactor Mr. 
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The fresh feed flow rate F 0 to the reactor and the reactor effluent flow rate F are 
now the unknown quantities, which must be computed from the model equations. 
The dimensionless time in the model governing the dynamic behavior of the system 
(obtained from Equations (3.2)-(3.11)) is now defined with respect to the recycle stream 
rate entering the reactor. At steady-state, these equations can be rearranged to yield 
the following single equation in 2 


(*2(z,P,n) = B(z-x e )(l-y e )-[y e -x e +p 0 Ti(y e -z)\ In 


(z-x e )( 1 - y e ) 


= 0 (3.32) 


[ Da 0 Tiz{y e - z) \ 

ti in the above equation represents the quantity M R jL. The parameter set p represents 
the parameters (B, fio, Dao,x e ,y e ). The solution to the set of equations 


G 2 = G 2z = <7 222 = G 2ti = 0 (3.33) 

was computed numerically with z,ti,B and 0o as the unknowns and with several sets of 
values for Dao, x e and y e . It was observed that the solution to the pitchfork singularity 
lies in the feasible region of the unknowns. From the computations, it appears that 
there is no solution to the winged cusp singularity, i.e., there is no solution to the set 
of equations 

G 2 = G 2z = G 2zz = G 2t = G 2zt , = 0 (3.34) 

Therefore, for this set of specified variables, it appears the system has a pitchfork 
singularity at its most singular point. 

The hysteresis and isola varieties were computed numerically. A typical projection 
of the varieties in the ( B,fio ) plane is shown in Fig. 3.8. The bifurcation diagrams in 
the different regions are shown in Fig. 3.9. Atleast one steady-state exists for all values 
of ti as is to be expected since now the normal form would be a cubic. The multiplicity 
features of the steady-states for this system are again similar to that of a CSTR[23]. 

The generic guideline proposed by Luyben to avoid the snowball effect[29], is that 
for isothermal recycle systems, the recycle rate should be a specified variable. We find 
that this need not be true in the case of non-isothermal recycle systems, though for this 
choice of the specified variables, feasibility boundaries do not influence the multiplicity 
features of the system. This is because fixing the recycle flow rate, constrains the 
conversion in the reactor to always lie within the bounds imposed by the vapor and 
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Figure 3.8: Global classification of the steady-state bifurcation diagrams for {T F ,P) 
and {M r , L ) flash (Equation (3.32)) in the B - /?<, plane for x e = 0.95, y e = 0.1, 
Da 0 = 0.10. HV: Hysteresis variety; IV: Isola variety. 

*1 





Figure 3.9: Schematic bifurcation diagrams of the steady-state (T F , P ) and {Mr, L) 
flash (Equation (3.32)) describing the dependence of composition in the reactor on t/. 
Letters indicate regions in Fig. 3.8. 
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liquid stream compositions exiting the separator. The reactor conversion z cannot 
equal y e , since this would mean that the recycle rate goes to zero. Neither can the 
reactor conversion z equal x e , since this would result in the fresh feed rate going to 
zero. As explained for the case (( Mr,F ) separation), this results in a situation where 
we do not have any product stream exiting the system. As a result, we have the system 
operating in batch mode, leading to an increase in the conversion and hence a decrease 
in z. Such a situation again arises only because we have already fixed the operating 
conditions (Tp, P) in the separator. 

However, in the case of z — >■ y e , i.e., at high conversions, the amount of recycle 
going into the liquid fraction decreases. Hence for a small decrease in z, the amount 
of fresh feed added to maintain the recycle rate at its specified value, increases. This 
results in a decrease in the residence time in the reactor. In the case of non-isothermal 
reactions, this decrease in the residence time in the reactor can result in the appearance 
of run-away like conditions. Another effect of such an interaction is in the appearance 
of the snowball effect, since a small decrease in z results in a very large increase in the 
fresh feed flow rate to the reactor. 

3.3.4 Isothermal and Isobaric Separation with (Fq, F) Fixed 

We now consider the case when both the fresh feed flow rate Fq and the reactor ef- 
fluent flow rate are the specified variables. Since we assume isothermal and isobaric 
separation, the material balance and equilibrium equations governing the flash (Equa- 
tions (3.5), (3.6) and (3.8)-(3.11)) determine the composition z of the reactor effluent 
stream. Hence this composition can be treated as a bifurcation parameter. The molar 
holdup Mr, the temperature T in the reactor, and the recycle flow rate L are now 
unknown. The dimensionless time in the model governing the dynamic behavior of 
the system (obtained from Equations (3.2)-(3.11)) is now defined with respect to the 
(unknown) actual residence time in the reactor. At steady-state, these equations can 
be rearranged to yield the following single equation in T 

(?3 ( T , p, z) = B(1 — yejDciQZ^Z X e ) T ^Z)QoZ(y e %e) “b A) (I He) (2 x e )e j = 0 

(3.35) 

The reactor conversion z is now considered to be the bifurcation parameter. This 
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implies that we operate the reactor at a specified conversion. The parameter set p 
represents the parameters (B,(3o,Da 0 ,x e ,y e ). The solution to the set of equations 

C?3 = Gzt = Gztt = G 3z = 0 (3.36) 

is given by 

X e (l - y e ) R _ 2(z e -Ve) R _ gMge ~ i/e)^ 

T = 2 ’ 2 = T^T’ b -^T^ 7' i-!/ 

This solution to the pitchfork singularity lies in the feasible region. Further, it can be 
shown that there is no solution to the winged cusp singularity, i.e., there is no solution 
to the set of equations 

Gz = Gzt — Gztt — Gzt — Gztz — 0 (3.37) 

Therefore, for this set of specified variables, the system has a pitchfork singularity at 
its most singular point. 

A typical projection of the hysteresis and isola varieties and the boundary sets in the 
( B,/3q ) plane is shown in Fig. 3.10. The bifurcation diagrams in the different regions 
are shown in Fig. 3.11. It can be seen from the figure that atleast one steady-state 
exists for all feasible values of the bifurcation parameter z. The multiplicity features 
of the steady states for this system are essentially similar to that of a CSTR[23]. As 
mentioned earlier, for this set of the specified variables, the conversion z and all the flow 
rates to and from the reactor are determined by the material balance and equilibrium 
relations for the flash. Since we are considering isothermal and isobaric separation, the 
bounds on the conversion z still exist. Feasibility boundaries in this case arise due to 
the following bounds on the bifurcation parameter 

y e < z <x e 

The two vertical lines in the bifurcation diagrams in Fig. 3.11 represent the feasibility 
boundaries corresponding to y e and x e (x e > y e ). Specifying the rates of the fresh feed 
stream to the reactor and the reactor effluent stream, implies that z cannot equal x e . 
Therefore, only the feasibility boundary corresponding to z = y e exists for this case. 
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Figure 3.10: Global classification of the steady-state bifurcation diagrams for ( Tf,P ) 
and (Fo, F) flash (Equation (3.35)) in the B - j3 0 plane for x e = 0.95, y e = 0.1, 
Da 0 = 0.10. BLS: Boundary limit set; HV: Hysteresis variety; IV: Isola variety. 
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Figure 3.11: Schematic bifurcation diagrams of the steady-state for (7>, P) and (F 0 ,F) 
flash (Equation (3.35)) describing the dependence of temperature in the reactor on 
z. Letters indicate regions in Fig. 3.10. The vertical lines represent the feasibility 
boundaries z = x e and z = y e respectively ( x e > y e ). 
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This situation corresponds to the coupled system behaving like a system in series, i.e., 
without recycle. The boundary limit set (BLS) for z = y e \s shown in Fig. 3.10. The 
effect of this boundary oh the multiplicity features of the coupled system is shown in 
Fig. 3.11. The part of the bifurcation diagram to the left of the z — y e line (the vertical 
line to the left) is lies in the infeasible region and hence is of no physical significance. 
This corresponds to a negative value of one of the flow rates. At the limit of z x e , 
the value of T approaches zero. This corresponds to a situation where we have very 
low residence times in the reactor and a very low amount of fresh feed entering the 
reactor. 

Practical considerations for a smooth operation of the coupled system mean that 
the system has to be operated at low values of B and high values of /3q. Selection 
of these parameters such that the bifurcation diagrams are of the type shown in Fig. 
3.11(a), (b) or (g) ensures that the feasibility boundary does not influence operation 
of the system at the steady-state. 

3.3.5 ( Tp,L ) Separation with (Mr, Fo) Fixed 

We finally consider the case when the liquid rate L leaving the flash is the specified 
variable in addition to the fresh feed flow rate F 0 and the molar holdup in the reactor. 
The reactor and mixer holdups are assumed to be constant. Therefore, the reactor 
effluent stream rate is now a known quantity. This enables us to define the dimension- 
less time in the model governing the dynamic behavior of the system (obtained from 
Equations (3.2)-(3.11)) with respect to the actual residence time in the reactor. The 
pressure P in the separator is now an unknown quantity which must be computed by 
solving the set of model equations simultaneously. At steady-state, these equations can 
be rearranged to yield the following single equation in x e 

(1 - x e )(a + 1) 

Da 0 Tf[x e (a + 1) + ctxK 7 - 1)] 

(3.38) 

In the above equation, 7 is defined as P{ at [P$ at . 77 represents the actual residence 
time in the reactor. The parameter set represents the parameters (B,/3 0 , Da 0 ,7,a). 
The solution to the set of equations 

C?4 = Gix e ~ G 4x e x e = - : G^ T j = 0 


G A (x e ,p,T f ) = B(l—x e )— [l+x e (7— 1)] [l+a+PoTf\ In 


(3.39) 
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Figure 3.12: Global classification of the steady-state bifurcation diagrams for (T F ,L) 
flash (Equation (3.38)) in the (a, ft) plane. 

was computed numerically with x e , tj , a and 7 as the unknowns and with several 
sets of values for B, Da 0 , and ft. It was observed that the solution to the pitchfork 
singularity lies in the feasible region of the unknowns. It can be shown that there is no 
solution to the winged cusp singularity i.e., there is no solution to the set of equations 

G 4 — G^Xe — GfocXc = Giy = Gfa'y = 0 (3.40) 

From the computations, it appears that for this set of specified variables ((7>, L) flash) 
the system has a pitchfork singularity at its most singular point. 

The hysteresis and isola varieties were computed numerically. A typical projection 
of the varieties in the (a, ft) plane is shown in Fig. 3.12. The bifurcation diagrams 
in the different regions are shown in Fig. 3.13. Atleast one steady-state exists for all 
values of 77. This is to be expected since now the normal form would be a cubic. The 

multiplicity features of the steady-states for this system are again similar to that of a 
CSTR[23j. 

In the present mode of operation, the choice of the specified variables is such that 
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Figure 3.13: Schematic bifurcation diagrams of the steady-state for (TV, L) flash (Equa- 
tion (3.38)) describing the dependence of temperature in the reactor on 77. Letters 
indicate regions in Fig. 3.12. 
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there are no pre-determined composition constraints acting on the system, since the 
pressure in the separator is an unknown quantity. We can visualise the system as 
having to evolve to a state, which satisfies the constraint imposed by fixing the recycle 
flow rate to the reactor. For example, the problem can be consideied to be one of 
determining the composition of the feed stream to the separator and the pressure in 
the separator, such that the flow rate of the liquid fraction exiting the flash equals 
its specified value. The absence of feasibility boundaries on the reactor conversion 
implies that the coupled system can regulate itself in this manner for all combinations 
of the operating parameters. For this reason we have depicted the dependence of the 
temperature in the reactor in the bifurcation diagrams. 

Investigating the different modes of operation of the system enables us to explain 
the reason for the non-existence of steady states observed in the (7>, P) and ( Mr , F 0 ) 
flash operation. The fundamental difference between this and the other modes can be 
best understood by considering the overall system performance. 

For the case of isothermal and isobaric separation when (Mr, F 0 ) are the specified 
variables, it was seen that the interactions between the two units (through recycle) is 
not allowed to affect the performance of the coupled system. This is because speci- 
fying the above four variables fixes the product rate and the product quality of the 
system. Therefore, this choice of combination of the specified variables constrains the 
performance of the coupled system. Non-existence of steady-states for certain values of 
parameters can be viewed as a failure of the coupled system to satisfy the constraints. 

For all other combinations of the specified variables, either the product flow rate 
or its quality is constrained. Therefore, for all these combinations of the specified 
variables, the interactions between the two units (through recycle) can affect the per- 
formance of the coupled system, through the unconstrained variable. The system is 
thus free to evolve, and does evolve, to a state such that all the constraints are satisfied. 
Therefore, we do not have any regions in parameter space where both the feasible and 
infeasible steady-states of the single steady-state equation are non-existent. In these 
cases, the nature of the imposed constraints determines the presence or absence of the 
feasibility boundaries. 

Interpretation of the results presented in this section can also be undertaken from 
the mathematical point of view. Absence of steady states for the case discussed in sub- 
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section 3.3.1 arises because the equation (3.14) is contact equivalent to a quadratic. All 
the other cases considered (sections 3. 3. 2-3. 3. 5) have a pitchfork singularity. Singular- 
ity theory predicts that the normal form for each of these cases is cubic in nature. This 
assures us of the presence of atleast one steady-state for all operating conditions. How- 
ever, the presence of bounds on the conversion in the reactor can result in the existence 
of regions in parameter space where we do not have any feasible steady states. 

In Figs. 3.7 and 3.9, for the cases of (Mr, F ) flash (subsection 3.3.2) and (Mr, L) 
•flash (subsection 3.3.3), we see that for very low residence times in the reactor, the 
steady-state conversion z approaches x e . Having z — > x e in these cases, corresponds to 
a situation where we have a very low amount of fresh feed entering the reactor. 

The analysis carried out in this section shows that it is possible to cross the z = y e 
boundary for some sets of specifications. Operation near this boundary this can lead 
to unsafe conditions, since this boundary would typically correspond to a state of high 
conversion and high temperature in the reactor. 

In the above discussion, we have concerned ourselves only with analysing the steady- 
state behavior of the system. The dynamic behavior has not been considered at all. 
The nature of the dependence of the steady-state composition and the temperature in 
the reactor (represented in the bifurcation diagrams) on an independent variable can 
be considered to be obtained by slowly varying the bifurcation parameter. Considering 
such a quasi-static variation of the control parameter together with the arguments 
presented earlier in this section enables us to obtain useful information about the 
dynamic behavior of the system. 

During its approach to the steady-state, i.e., during the transient evolution, the 
system can evolve towards and cross the feasibility boundaries represented by x e or y e . 
As discussed earlier in this section, such an evolution of the system can result in the 
failure of operation of the coupled system, if appropriate corrective action is not taken. 

3.4 Conclusions 

In this work we have investigated the behavior of a coupled reactor-separator system. 
The steady-state multiplicity features were studied using singularity theory. It was 
shown that the choice of the specified variables drastically changes the multiplicity 
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features of the system. It was shown that for a paiticular choice of the specified 
variables (( M r ,F 0 ,T f ,P ) operation), the maximum number of feasible steady-states 
of the system is two. For this choice of the specified variables, theie exist ianges of 
parameter values for which a steady-state of the coupled system can be non-existent. In 
this case, the non-existence of the steady-states arises because the performance of the 
coupled system (in terms of production rate and product composition) is completely 
determined by the VLE considerations and the overall mole balance of the coupled 
system. The performance of the coupled system is not allowed to be influenced by the 
dynamics and interaction between the individual units. In contrast, the multiplicity 
features with the other sets of specifications, are essentially similar to that of the CSTR, 
which has a maximum of three steady-states. The only difference is in the presence of 
additional multiplicity features resulting from the presence of feasibility boundaries. 

An important feature of this system is that the multiplicity features are influenced 
by feasibility constraints generated by the VLE considerations in the flash. This is 
true only for isothermal and isobaric separation. When either the temperature or the 
pressure in the flash are not specified, the VLE relations do not impose any constraints 
on the system. Therefore, in this case, eventhough the some of the flow rates are 
constrained, the vapor composition can vary since the pressure in the flash is not fixed. 
In this case the system can regulate itself, such that the flow constraints are satisfied. 
As a result, we observe that the feasibility boundaries on the conversion do not exist. 

It is a common practice in industry to operate a distillation column by maintaining 
the compositions of the top and bottom streams at constant values. Operating a reactor 
with constant fresh feed flow rate is also very common. This control configuration, 
however, may not be the best option for a recycle system. Our results indicate that 
for a simple recycle system consisting of a reactor and a separator, such an operation 
can lead to severe problems in the operation of the system. 



Chapter 4 

A Unified Framework for the 
Analysis of Synchronization of 
Chaotic Systems 


Nomenclature 


Pi, Pi, P3 

Inputs in the Rossler system. 

r 

Relative order of the system. 

t 

Time. 

u 

Manipulated input. 

x m 

States of the reference system. 

Xp 

States of the process system. 

y 

Process output. 

Vd 

Desired or reference output. 

z 

Co-ordinates in the transformed co-ordinate system. 


4.1 Introduction 

The temporal evolution of systems which exhibit chaotic behavior is characterized by 
a sensitive dependence on initial conditions. This implies that the trajectories along 
which two identical systems starting from nearby initial states evolve, diverge from 
each other. Mathematically, such an evolution is characterized by a positive value of 
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Drive system 


Figure 4.1: Schematic diagram of Pecora-Carroll synchronization. 

one or more Lyapunov exponents. The Lyapunov exponents, which characterize the 
nature of a trajectory along which a system evolves, are analogous to eigenvalues which 
determine the stability of a steady-state of a system. Positive (negative) values of the 
Lyapunov exponents, are a measure of the rate at which two nearby states diverge from 
(approach) each other, as the system evolves. 

Synchronization refers to the ability of all the states of two identical systems, start- 
ing from nearby initial states, to evolve identically with time. It would hence appear 
that two chaotic systems cannot be made to synchronize with each other. However, 
Pecora and Carroll[37, 38] showed that two chaotic systems could be synchronized with 
one another. They achieved synchronization between two chaotic systems by linking 
them with a common signal. A schematic representation of their technique is shown 
in Fig. 4.1. They considered a three-dimensional system represented by the state vari- 
ables X , Y and Z. The state variables of the first and second systems are denoted by 
the subscripts 1 and 2 respectively. System 1, characterized by X x , Yj and is called 
the drive system. System 2, characterized by X 2 , Y 2 and Z 2 is called the response sys- 
tem. Pecora and Carroll studied the evolution of the response system , represented by 
the states Y 2 and Z 2 , when X 2 was equal to X u i.e., when the Y 2 , Z 2 system was driven 
by the X x signal. Mathematically, this can be represented by the following system of 
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Figure 4.2: Schematic diagram for detecting the occurrence of generalized. synchroniza- 
tion in a given system. 

equations 

X x = F 1 (X 1 ,Y 1 ,Z 1 ) 

Yi = F 2 (X 1 ,Y 1 ,Z 1 ) Y 2 = F 2 (X u Y 2 ,Z 2 ) 

Zi = F 3 (X 1 ,Y 1 ,Z 1 ) Z 2 = F 3 (X u Y 2 ,Z 2 ) 

Since the driving signal Xy is chaotic, the evolution of the response system is also 
chaotic. The drive and the response systems are said to have synchronized, if the 
difference between the corresponding states of the two systems,- i.e., (Yi—Y 2 ) and (Zy — 
Z 2 ), goes to zero as t — * oo. Pecora and Carroll showed that synchronization of the drive 
and response systems was possible, if the response system has only negative Lyapunov 
exponents, As mentioned in the first chapter, the importance of synchronization stems 
from its possible applications in several disciplines of science. 

Recently, two different types of synchronization have been reported in the literature. 
Rulkov et al. [39] and Kocarev and Parlitz[40] have studied a type of synchronization, 
which they refer to as generalized synchronization (GS). The driving signal (Fig. 4.2) is 
now a function (possibly nonlinear) of the states of the drive system, i.e., it is achieved 
by substituting the driving variable (X 2 ) with the driving signal in the response system. 
Using such a driving signal generally results in a nonlinear relationship between the 
states of the driving system and the corresponding states of the response system. A 
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Figure 4.3: Schematic diagram of the feedback control strategy. 

simple way of detecting whether or not synchronization is achieved involves construct- 
ing a replica of the response system (Fig. 4.2). The response system is said to have the 
ability to synchronize if the errors (V 2 - Y 3 ) and (Z 2 — Z 3 ) go to zero as t — > 00 . Ko- 
carev and Parlitz and their co-workers[41, 42] have reported a powerful new technique 
for achieving synchronization in chaotic systems. Their procedure, an active-passive 
decomposition (APD) technique, involves creating a new system by defining a new 
variable as a function of the drive system variables. The new variable is used as the 
driving signal to drive the modified response system. This is discussed in greater detail 
in the third section of this chapter. 

In this chapter we present a general framework for the analysis of the above tech- 
niques for synchronizing chaotic systems. We show that the above techniques are a 
different interpretation of the use of feedback, to achieve perfect control of a process 
along a desired (chaotic) trajectory. By perfect control, we mean that the controlled 
output tracks the desired output exactly for all time t > 0 . We show that the three 
types of synchronization mentioned above are different cases of perfect control of a sin- 
gle input-single output (SISO) system. The connection between synchronization and 
perfect control is illustrated on a modified form of the Rossler system. 
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4.2 Theoretical Analysis 

Figure 4.3 shows a schematic diagram of the feedback control problem. System 1 
represents a model whose output y<i generates the desired (reference) behavior of the 
process. System 2 represents the actual process with an output y. We assume that 
the model represents the process exactly. The difference between the evolution of 
the two systems arises due to the differences in the initial conditions of the state 
variables x m and x p . The process output y cannot track the desired output yd in the 
absence of control, because of the chaotic behavior of the system. It is desired to match 
(synchronize) the process output y along a desired trajectory yd for t > 0 using output 
feedback control. This objective is different from the conventional control goal, where 
it is desired that synchronization occurs as time t -¥ oo. 

Analysis and design of a controller for achieving perfect control in a nonlinear 
system requires the use of results from nonlinear systems and nonlinear control theories. 
These theories for nonlinear systems which are linear in the input or manipulated 
variable are well developed. They are an extension to nonlinear systems, of the Laplace 
transform method for linear systems. Recent developments in differential-geometry 
based nonlinear control theory [46, 62, 63], have made it possible to devise control 
schemes for the perfect control of nonlinear systems. We next present the formulation 
of the control problem and its solution. Complete details are available in[62, 63]. 

We consider the desired output and the process to be represented by the following 


n-dimensional single input-single output (SISO) models 

Xm = f(x m ) + g(x m )u 0 (4.1) 

yd — h(x m ) (4-2) 

Xp = f(xp) + g(x p )u (4.3) 

V = ft(x P ) (4.4) 


The overdot indicates differentiation with respect to time. x m ,x p € R n are the ref- 
erence and process states, u 0 and u are scalars representing the nominal and the ma- 
nipulated inputs respectively. For the purpose of showing the equivalence between 
synchronization and perfect control, it is sufficient to assume that only one element in 
the vector g(x p ) is non-zero. This implies that the input u, directly affects only one of 
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the states of the system, y is the scalar measured output, and yd denotes the scalar de- 
sired output to be tracked. Our objective is to apply suitable parametric perturbations 
on u, such that the process output y tracks the desired output yd exactly, i.e., 

Vd~y = 0 for ' t > 0 

Differentiating (4.4) repeatedly, the relationship between the output y and the state 
variable x p can be expressed as 

y = M* P ) 

* = L ' h M 


(4.5) 


- L r f l h(x p ) 

S? = qh(x v ) + L s Li- l h(x p )u 

The index r is a measure of how directly the input u affects the output y. It is called the 
relative order of the system and represents the number of differentiation operations 
which must be performed on the output, in order to obtain an explicit dependence of 
the output on the input u. The relative order is defined as the smallest integer r such 
that 


LgLf ^(xp) ^ 0 

where Lf/i(x p ) is the Lie derivative of a function h(x p ) with respect to f (x p ) and is 
defined as 


Higher order Lie derivatives with the same vector argument are defined recursively as 

Tf 1 = TfLjr 1 
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The relationship between the 
be expressed as 


manipulated input u and the process output y, can then 


- L r f h(x p ) 

L s L r f l h(x p ) 


(4.6) 


where y r is the r th order derivative of the process output. For the problem under 
consideration, the measured and the desired outputs are specified by the designer. In 
order to obtain the control law which achieves the objective of perfect control, we set 
the output tracking error e = y — yd equal to zero. On differentiating this r times, we 
obtain the following equation relating the input u to the output tracking error e. 


e r = y r - y r d = 0 

The control law which achieves the objective of perfect control is obtained by solving 
the above equation for u. Equivalently, this is given by replacing y T in (4.6) by y r d . It 
follows that, by implementing the control law, the desired trajectory can be tracked 
perfectly, i.e., y(t) = y d (t ) for all t > 0. 

The tracking error equation above is of order r. In most cases r < n, and the 
tracking error equation accounts only for a part of the dynamics of the closed-loop 
system. The effect of this is to render part of the process dynamics unobservable 
from an input-output point of view. This part of the system dynamics is called the 
interned dynamics of the system, because the input-output relationship provides no 
information about the behavior of the states representing the internal dynamics of the 
system. The control is successful only if the internal dynamics are stable. Stability of 
the system representing the internal dynamics can be obtained by defining the zero 
dynamics of the system. The zero dynamics are an intrinsic feature of nonlinear 
systems and are defined as the internal dynamics, when the process output is kept at 
its desired value by the input. 

Information about internal stability of the process (i.e., whether all the process 
states of that part of the system representing the internal dynamics, also approach 
their desired trajectories), can be obtained from a stability analysis of the zero dy- 
namics of the process. Systems with asymptotically stable zero dynamics are said to 
be minim um -phase systems. Next we outline the method for obtaining the zero 
dynamics of a system. 
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Generalization of the analysis, for obtaining the zero dynamics is possible by work- 
ing in transformed co-ordinates. Since Lf h(x p ), k = 0, • • • , r — 1 are linearly indepen- 
dent functions, we can choose these to be the first r elements in defining a transformed 
co-ordinate system z, i.e., 

z k = L k f l h(x p ) /c = 1, 2, ■ • • , r 

It is also possible to define (n - r ) additional co-ordinates, Zk = 0fc(x p ),fc = (r + 
1), ■ • • , n, such that their time derivatives are independent of u[46, 64]. These additional 
co-ordinates can be obtained from the (n — 1) solutions w(x) of the partial differential 
equation 

ft (*p) gr + ft(xp) £ + ••: + «»(*p) <£; = 0 

The (n — 1) solutions to the above partial differential equation are obtained by defining 

%p,i+l = Cl) ' ' ' ) Cn— l)) i — 1) " ‘ ’ ) 11 1 (4.7) 

to be the solution of 

= **(«)= G. i = V",n 

Then the equations in (4.7) above, are solvable for (i, • • • , C„_i, i.e., 

C i = &(x p) 

Further, the scalar fields Zk = <^fc(x p ), k = 1, 2, ■ • • , (n — r) thus obtained, are linearly 
independent and satisfy the above partial differential equation. 

In the z co-ordinates, the SISO nonlinear system can be represented in its normal 
form 

Zl = *r+l = ?r+l(z) 

i r = b{z) + a{z)u z n = q n ( z ) 

y = zi 


where 


a(z) = L g Lf ^[r 1 ^)], 5(z) = L r { h[<f>-\ z)] and 

9fc( z ) = ^z)] A: = 1, • • • , (ji — f) 
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A minimal-order realization of the inverse (MORI) of the system is obtained by re- 
placing u in the normal form above by (4.6), after a suitable change of co-ordinates. 
The advantage of using the normal form representation of nonlinear SISO system is 
that its inverse is effectively of dimension (n — r). The zero dynamics of a system are 
defined as the dynamics of a MORI. As mentioned earlier, they represent the internal 
system dynamics when the process output is constrained to be the desired output. The 
dynamical system representing the internal dynamics of the system is given by the last 
(n — r ) equations of the normal form. 

Before showing the equivalence of synchronization and control with examples, we 
explain how the three types of synchronization, result from consideration of the control 
problem. Examples are presented in the next section to show that synchronization is 
another interpretation of perfect control when r = 1. 

The synchronization technique proposed by Pecora and Carroll[37] results from the 
case when both the desired and measured outputs are states (yd = x mi and y = x pi , for 
an i € (l,n)). GS[39, 40] results from considering the desired output as a function of 
states yd = hi(x m ) and y = h 2 (x p ) = x pi , 1 < i < n. For the case of synchronization 
using the APD technique, we rewrite one of the system variables (after neglecting the 
subscripts) as 

ii = fi(x) + gi(x)u + e(x) - e(x) 

where e(x) is an arbitrary function of the states x. We can define a new variable x n+ i, 
as the complete (or part of) function /;(x) + gi(x)u -1- e(x). Further, we can represent 
the evolution of this new variable by an arbitrary equation, containing a new input 
u new . The model and process can now be represented by equations of the form 

x = f(x) + g(x)x n+ i - e(x) 

Xn+i = /i(x) +gi(x)u + e(x) 

Xn +1 = d(x) 4" u new 

Now consider the case where the desired and process output are given by ya = x m , n +i 
and y = x Pjn+ i, and the manipulated input is u new . Perfect control, for this case, is 
the control analog of synchronization using APD. From the point of view of synchro- 
nization, the exact functional form of d(x) is not required since the desired trajectory 
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Vd - x m ,n+u is directly obtained from a knowledge of the model states x m . Several 
examples of this case can be found in Table 1 of Pailitz et, <zi.[42]. 

4.3 Results and Discussion 

In this section we demonstrate the connection between perfect control and synchroniza- 
tion with the help of an example. The Rossler system is one of the extensively studied 
models in the nonlinear dynamics literature. Further, different aspects of chaotic syn- 
chronization of this system have been studied[37, 38, 42]. For these reasons, we have 
considered this system in this work. Our objective will be design the controller such 
that the process output tracks the the desired output exactly. This enables us to 
show the equivalence between the three different types of synchronization and perfect 
control of a process of relative order one. As mentioned in the previous section, we 
consider that the reference system models the process system exactly. We represent 
the reference system (system 1) and the process (system 2) as 

Xi = -Vi-Zi (4.8) A 2 = -Y 2 -Z 2 + Pi (4.11) 

Fi = 4- aYi (4-9) Y 2 — X 2 -1- ( 1 Y 2 + p 2 (4.12) 

Zi = b + Z t {Xx-c) (4.10) Z 2 = b + Z 2 (X 2 - c) + p, (4.13) 

The system parameters a, b and c are taken to be the same for both the systems. We 
have chosen a — b = 0.20 and c = 9 for our studies. The system exhibits chaotic 
behavior for these values of the parameters. We have introduced p \ , p 2 and pz as 
additively occurring auxiliary inputs in the process system. These are used as the 
manipulated inputs for achieving the control objective. We also assume that only one 
of the auxiliary inputs (pi,P 2 ,Pz) is non-zero. This is the manipulated input. 

Consider the case where y = X 2 , u = Pl (with p 2 = p z = 0) and y d = X x . The 
manipulated input p\ occurs in the governing differential equation of the controlled 
output X 2 and hence r = 1 . For perfect control, we require that y — y d for t > 0 . This 

means that at t = 0 , we require X 2 — X x . The control law for tracking the desired 
trajectory y d is 

u = Pi=y d - (~Y 2 - Z 2 ) ( 4 . 14 ) 
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Substituting this back in the equations of the process (4.11)-(4.13), we get 


X 2 

= ijd 

( 4 . 15 ) 

y 2 

= yd + aY 2 

( 4 . 16 ) 

^2 

— b + Z 2 (yd — c ) 

( 4 . 17 ) 


Applying the control action computed using the control law (4.14) ensures that the 
desired output can be tracked exactly, i.e., y = y d and X 2 = X x , for t > 0. However, 
this provides no information about the behavior of the remaining system states, Y 2 and 
Z 2 . Information about the behavior of these states can be obtained by studying the 
stability of the internal dynamics of the process system. We now address the issue of 
internal stability of the process. 

The system (4.11)-(4.13) is in its normal form. The co-ordinates z of the system in 
its normal form are the process system states. Hence the states Y 2 and Z 2 constitute 
the states of the internal dynamics of this system. The dynamical system representing 
the internal dynamics is given by (4.12) and (4.13). The zero dynamics are obtained by 
replacing y = X 2 in these equations by yd (Equations (4.16) and (4.17)). The stability 
of the zero dynamics can be obtained by carrying out a stability analysis of the zero 
dynamics of the system. Since the desired trajectory is chaotic, the necessary condition 
for asymptotic stability of the zero dynamics (Equations (4.16) and (4.17)) is that its 
Lyapunov exponents must be negative[37]. These were calculated to be +0.20 and 
—8.87. This indicates that the zero dynamics is unstable (the system is non-minimum 
phase), and hence the process is not internally stable, i.e., the variables Y 2 and Z 2 do 
not approach their respective trajectories. This can also be seen by taking the error 
between the states Y\ and Y 2 . Defining e 2 to be the error (Yi — Y 2 ), we have 


e 2 — a(Yi - Y 2 ) = ae 2 


This system , representing a part of the internal dynamics of the process system, is 
unstable. The Lyapunov exponent of this part of the internal dynamics is a, i.e., 0.20. 
The growth of the error e 2 with time is shown in Fig. 4.4. 

The set of equations (4.8)-(4.10), (4.16) and (4.17), is identical to the case of X 
drive in the approach of Pecora and Carroll. This is because the control parameter p x 
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Figure 4.4: Evolution of the state error e 2 = Y l - Y 2 with 
Vd — Xi and u = p 1 . 
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Figure 4.5: Evolution of the state errors — (X\ — X 2 ) and e 3 = [Z\ — Z 2 ) with time 
for the case of y — Y 2 , yd = Y\ and u = p 2 . 
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Figure 4.6: Evolution of the error e x = (X x - X 2 ) and e 3 = (Z x - Z 2 ) with time for the 
case of y = Y 2 , yd = Yi and u = p 2 . Parameter mis-match of ±2% between the model 
and the process was assumed. 


important to study the dynamics of the response system, by considering the effect of 
parameter mis-match. We illustrate the results for the case of y = Y 2 , yd = Y 2 and 
u — p 2 {Y drive). Figure 4.6 shows the evolution of the errors e\ — (Xi — X 2 ) and 
e 3 = [Z\ — Z 2 ) when there is a +2% and —2% mis-match between all the parameters of 
the model and the process. The reduction in the quality of synchronization is clearly 
seen from the figure. Some amount of off-set appears to exist for the e x component of 
the error. 

Next we show the equivalence between perfect control and GS[40]. We consider the 
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following control problem 

y, $ = X\ + Y\ Zi, y — X 2 and u — p\ 

This corresponds to synchronization with the driving vaiiable being a function of the 
states of the drive system. For this case, r = 1 and the system (4. 11 )-(*!. 1.3) is in 
its normal form. As mentioned earlier (discussion following Fig. 4.4), control (syn- 
chronization) cannot be obtained with this choice of the output (drive), since the zero 
dynamics (Equations (4.16) and (4.17)) are unstable. For the Rossler system, it is 
possible to achieve GS using the same desired output but by considering y = Y 2 and 
u=p 2 . 

For an example of synchronization using the APD technique, we begin by defining a 
new variable x n+ i — s = 1.2 Y. Using this, the Y component of the model (4.8)-(4.13) 
can be rewritten as 


Y 1 = X l -Y l + s l and Y 2 = X 2 -Y 2 + s 2 (4.18) 

This new variable is selected such that the response system consisting of the states X 2) 
Y 2 and Z 2 , is asymptotically stable. We now assume that the evolution of this new 
variable is given by an equation of the form 


in+i = s= d(x) + u new (4.19) 

where d(x) is an arbitrary function of the states x. The dynamical system representing 
the reference and the process systems obtained by incorporating this new variable can 
then be obtained using the above equations for Y\ and Y 2 . The perfect control objective 
V — s 2 i Vd — s 1 = 1.2Ti with u — p 2 is the control analog of synchronization using the 
APD technique. The process system is in its normal form and the equations (4.11), 
(4.13) and the equation for the evolution of Y 2 given in (4.18), now constitute the 
dynamical system representing the internal dynamics of the modified system. 

We integrate the resulting system of equations by assuming s 2 — Si. This obviates 
the need to specify the function d(x). Figure 4.7 shows the evolution of the state 
errors with time. When the model represents the process exactly, the tracking error 
goes to zeio very rapidly. This is shown in Fig. 4.7(a). However the performance of the 
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designing a nonlinear controller for a system of relative oitlci gieatei than one, i.e., we 
consider perfect control of the process when the manipulated input, u, does not directly 
affect the controlled output y. We illustrate the method by considering A'-, to be the 

controlled output (y = W 2 and y d = Xi) and p 2 to be the manipulated input (u = p 2 ). 

An explicit relationship between the input and the output can be obtained after per- 
forming two differentiation operations on the output. Hence we have r = 2 for this 
case. Proceeding as described in the previous section, the control law is obtained as 

u — p 2 — -yd - Z 2 - xjd - o.Y 2 (4.20) 

Substituting the control law back into the equations describing the process (4.11)- 
(4.13), we get 

X 2 = y d - X\ (4.21) 

y 2 = -y d -Z 2 = Y 1 + Z l - Z 2 (4.22) 

Z 2 = b + Z 2 {y d ~c ) (4.23) 

Requiring perfect control on initiating the process setting A r 2 = A'j (?/ = y d ) at t = 0 
requires that X 2 = X x for t > 0. On rearranging, this yields: Y 2 — —y d - Z 2 . This 
implies that manipulating p 2 in order to make X 2 track X\, constrains Y 2 to evolve 
such that the above relation is satisfied. Hence it is only possible to achieve perfect 
control if at t = 0, Y 2 = —y d — Z 2 . The co-ordinates z of the process system in its 
normal form are related to the actual co-ordinates by the relations 

X 2 , z 2 — Y 2 Z 2 , z 2 — Z 2 

The normal form of the system is represented by the following system of equations 
z i = z 2, z 2 - b( z) + a(z)p 2 and i 3 = q( z) 

where a(z) = 1, 5(a) = + a (-z 2 - z 3 ), and q{z) = b + z 3 (z x - c) 

The internal stability of the process is now dependent on the asymptotic stability of the 
zero dynamics, given by (4.23). This system has a negative Lyapunov exponent equal 
to -8.87, indicating that sychronization between the process and the reference system 
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Pi 

P 2 

P 3 

*2 

0.20 

-8.87 

-8.87 

0.20 

y 2 

-8.87 

-0.024 

-8.8 


y 2 

0.20 

I 

0.10 

0.10 


Table 4.1: Lyapunov exponents of the zero dynamics for the process system for different 
control configurations of the Rossler system. 

can be obtained using feedback control. The evolution of the error e% = [Z\ — Z 2 ) with 
time is similar to that shown in Fig. 4.5. 

Lyapunov exponents of the zero dynamics for each combination of the input-output 
(output=state) are given in Table 4.1(a). These are analogous to the conditional 
Lyapunov exponents of[37, 38]. The diagonal cases correspond to synchronization 
using X , Y and Z drives in the approach of Pecora and Carroll. Non-diagonal cases 
for which we have only one Lyapunov exponent can be considered to be similar to 
synchronization with a two variable drive in the approach of Pecora and Carroll[37, 
38]. The absence of zero dynamics for the Y 2 -p$ and Z 2 -p 2 configurations indicates 
that the Lyapunov exponents are equal to zero, for these cases[63]. The implication 
of this is that in order to achieve perfect control with these combinations of the input 
and the output, we need to fix all the initial states of the process. 

4.4 Conclusions 

In conclusion, we have presented a nonlinear control technique, to achieve perfect 
control of a system along a desired trajectory. An important result of this, work is that 
it connects the idea of perfect control, with the synchronization phenomenon discussed 
in the physics literature. In this connection, we have presented results from the control 
literature, to show that chaotic synchronization can be analysed within the framework 
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of nonlinear control theory. 

Several methods have been proposed in the literature for the control of chaotic 
systems [6 5]. Very little attention has been given to providing a theoretical basis for 
most of these techniques. Also very little attention has been given to applying the 
results from conventional control theory to the problem at hand. The need therefore, 
appears to point towards applying results from control theory, to provide the theoretical 
framework for several of the methods proposed in the physics literature. This can also 
result in a common theoretical framework of techniques for control of chaotic systems. 
An advantage of this could be that a technique suitable for the problem at hand, could 
then be selected from the different options available. 



Chapter 5 

Synchronization and Control of 
N on- Minimum Phase Systems 


Nomenclature 

a Constant in Equations (5.17) and (5.30). 

k Dimensionless reaction rate constant for the autocatalytic reaction. 

t Time. 

u Manipulated input. 

w State of the exo-system. 

x m States of the reference system. 

x p States of the process system. 

y Process output. 

yd Desired or reference output. 

ft Dimensionless feed concentration of B in the CSTR. 

r Dimensionless characteristic time in the CSTR. 


5.1 Introduction 

In the previous chapter, we illustrated how results from nonlinear control theory can 
be employed to show the equivalence between synchronization of chaos and perfect 
control of a system along a desired trajectory. We also saw that not all systems can be 
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Figure 5.1: Schematic bifurcation diagrams for system showing input multiplicity. 

synchronized using the approach of Pecora and Carroll. An example of such a situation 
is when we use the X drive of the Rossler system[37, 38]. As we saw in the previous 
chapter, it is not possible to achieve with this drive since the zero dynamics of the 
system for this choice of the output (A) and the input pi is unstable. Systems with 

unstable zero dynamics are classified as non-minimum phase systems in the control 
literature. 

An important class of non-minimum phase systems of importance in chemical en- 
gineering are the systems exhibiting input multiplicity [43-45]. In these systems, it is 
possible to have more than one possible value of the control input, for a given value 
of the output. Figure 5.1 shows a schematic representation of two typical steady-state 
bifurcation diagrams for such systems. For the desired choice of the controlled output 
x (denoted by the horizontal line), we have two possible choices u t and u 2 of the ma- 
nipulated input u. Unmeasured disturbances acting on the plant can destabilize the 
controlled system from its original set-point (at Ul ) and cause the system to reach the 
other steady state (at u 2 ). Though the controlled output has the same value for each 
of these inputs, the values of the process states, corresponding to the value of the input 
u 2 will be different from the desired values. 

cent years, there have been major developments in the application of nonlinear 
techniques based on differential-geometry for the control of chemical engineering 
y , 62, 63], Experimental and computational studies have been carried out to 
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Figure 5.2: Schematic diagram of the feedback control strategy for regulating the 
output of a nonlinear system. 

develop and test a variety of techniques. A limitation of most of these techniques is 
that they are valid only for systems with stable zero dynamics, i.e., for minimum-phase 
systems. Hence they are of limited use for the control of systems which can exhibit 
input multiplicity. Recently Isidori and Byrnes[47] have developed a technique for 
the regulating the output of nonlinear systems. Their method is based on bifurcation 
theory and is suitable for the control of non-minimum phase systems. 

Their technique is aimed at controlling a nonlinear system (plant) in order to have 
its output track a reference signal generated by an external system (exo-system) [46, 
47]. The technique involves considering a composite system consisting of the plant and 
the exo-system (Fig. 5.2). The exo-system is selected such that its linearization has all 
eigenvalues on the imaginary axis. This choice ensures that the composite system has at 
least one eigenvalue on the imaginary axis. Such a choice of the exo-system guarantees 
the existence of a center manifold for the composite system. Isidori and Byrnes show 
that the regulator problem can be solved by selecting the feedback control law, such 
that the center manifold is rendered invariant, i.e., the trajectories which begin on this 
manifold, remain on it for all time. Once the center manifold of the composite system 
has been computed, a choice of the feedback control law is made such that this center 
manifold is rendered invariant. Selecting such a control law causes the output of the 
composite system to be zero at each point on the center manifold, thereby ensuring 
that the plant output is regulated. The only necessary condition for the solution of the 
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regulator problem is that the plant be exponentially stabilizable b} feedback. 

In the first part of this chapter, we describe how this technique can be used to 
formulate and solve the control problem. In the second pait of this chaptei, we apply 
the technique to achieve synchronization of the Rosslcr system along a desiied chaotic 
trajectory. We employ the technique to track a reference chaotic trajectory. Finally, 
we consider the cubic autocatalytic reaction in an isothermal CSTR[8]. We apply the 
method for the control of this reactor system which shows input multiplicity. 

5.2 Description of the Method 

We consider the SISO nonlinear process and the exo-system to be given by the following 
system of equations 


x = f(x) + g(x)u (5.1) 

w = s(w ) (5.2) 

The overdot indicates differentiation with respect to time. The first equation describes 
the plant with state x and input u. The second equation describes the scalar exo- 
system. f (x) and g(x) are vector functions. An equation relating the error e pt between 
the actual plant output y = h(x) and the reference trajectory = q(w) can be 
represented as 

e pe = h(x) + q(w) (5.3) 

It is desired to design the controller such that the error e pe -4 0 as t -4 oo. Finally we 
assume that (x, u) = (0, 0) is an equilibrium of the system. 

The problem of designing a state feedback controller for the plant involves deter- 
mining a control law of the form 


u = a(x,w) 


such that 


(5.4) 


• The closed-loop system has an asymptotically stable equilibrium at x = 0 when 
the exo-system is disconnected, i.e., when w = 0. 
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• For all initial conditions x(0),w(0) sufficiently close to the origin, e pe (t) -*■ 0 as 
t — y oo. 

The method proposed by Isidori and Byrnes [47] for constructing the solution to the 
regulator problem is based on the following assumptions. 

1. w = 0 is a stable equilibrium of the exo-system. This restricts the choice of the 
exo-system, to only those systems whose linear approximation at w = 0 has all 
eigenvalues on the imaginary axis. Selecting a stable exo-system ensures that 
w — > 0 as t — » oo. 

2. The pair f(x), g(x) has a stabilizable linear approximation at x = 0, i.e., if the 
linearization of the nonlinear system (5.1), (5.2) is given by 

x = Ax -I- Bn 
w = Sw 

where A = H| x =o, B = g(0) and S = There exists a matrix K such that all 
the eigenvalues of the matrix A + BK have negative real parts. 

Assumption 1 above implies that the composite system, consisting of the given system 
and the exo-system, has atleast one eigenvalue on the imaginary axis. Under these 
conditions, the composite system has a center manifold at (0, 0), given by the graph of 
the mapping x = 7r(w). Isidori and Byrnes[47] prove that the regulator problem can 
be solved, if and only if there exists a mapping x = n(w) with tt( 0) = 0 and 7r'(0) = 0, 
satisfying the conditions 

Q'j t 

— s(w ) = f(7r(w)) + g{7r(‘w))a(jr(w),w) (5.5) 

ow 

h(jr(w)) + q(w) = 0 (5.6) 

Computing the functions n(w) involves solving the above partial differential equation. 
It is not possible to obtain an exact solution in many cases. An approximation of the 
solution can be obtained by seeking the solution in the form of a power series[66] 


7 r(w) = ciw 2 -I- c^w 3 + c^w 4 -I 
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The polynomial coefficients are determined by equating terms of the same order on the 
left and right hand sides. 

Existence of the mapping x = n(w) is a direct consequence of the existence of a 
center manifold for the composite system. The feedback conti ol law is then determined 
such that trajectories starting on this center manifold remain in it for all time, i.e., this 
center manifold is rendered invariant. Thus our objective is to find the feedback law 
such that the output e pe of the composite system is zero at each point on the center 
manifold. 

The solution of the state-feedback regulator problem (SFRP) is given by a feedback 
law of the form 

a(x, w) = c(w) + K(x - 7r(iu)) (5.7) 

It can be shown that c(w) = a(7r(te), io)[47]. The second part of the control law is 
required, to bring the trajectories beginning away from the center manifold, on to it. 
Once this has been achieved, the first part of the feedback law ensures that the plant 
tracks the reference signal exactly. 

5.3 Results and Discussion 

We now illustrate the applicability of the method on the Rossler system. Figure 4.2 
shows a schematic diagram of the feedback control problem. System 1 is a model whose 
output yd generates the desired (reference) behavior of the process. System 2 represents 
the actual process having an output y. We assume that the model represents the process 
exactly. The only difference between the two systems, is in the initial conditions of 
the state variables x m and x p respectively. The process cannot follow the desired 
trajectory in the absence of feedback, because of the chaotic nature of evolution of the 
system. It is desired to match (synchronize) the evolution of the controlled process 
output along its desired trajectory as time i — > oo using state feedback. Additionally, 
we require stability of the closed-loop system when the control loop is implemented. 

At the moment, the control problem mentioned above is in a different form com- 
pared to the regulator problem (Fig. 5.2) discussed in the previous section. Therefore, 
in comparison with the previous chapter, we adopt a slightly different approach for reg- 
ulating the chaotic output of the Rossler system. We work with a model obtained by 
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taking the error between the corresponding process states x p and the reference syst ems 
states x m . The model in this form will be referred to as the plant. The variables in 
the model of the plant are the state errors e = x m — x p . In the absence of any control 
action, the error model has an unstable zero solution, consistent with the chaotic na- 
ture of evolution of the original system. Our objective will be to control this no nlin ear 
system (the plant), in order to have its output y p track a time- varying reference signal 
y e generated by the exo-system. The exo-system is selected such that it has a stable 
zero solution. When the control action is successful, the states e of the plant go to 
zero. This causes the actual states x p of the process to track their respective desired 
trajectories x m . 

The Rossler system is represented by the following set of equations 

X = -Y-Z (5.8) 

Y = X + aV (5.9) 

Z = b + Z(X-c) (5.10) 

Here a, b and c are constant parameters. Throughout this work, we assume a = b = 0.20 

and c = 9.0. As explained in the previous chapter, we represent the reference system 
(system 1) and the process system (system 2) with the following set of equations 


II 

1 

X 

1 

(5.11) 

X 2 

— — Y 2 — Z 2 -f- u 

(5.14) 

= Xx+aTi 

(5.12) 

y 2 

= X 2 + aY 2 

(5.15) 

= b + Z 1 (X 1 - 

-c) (5.13) 

z 2 

= b + Z 2 {X 2 -c) 

(5.16) 


We illustrate the case where the controlled output is a function of a single state. We 
consider X 2 to be the process output and u represents the manipulated input. As 
explained earlier in this section, implementation of the technique described in the 
previous section is facilitated by recasting the equations (5.11)-(5.16) in error form. 
These result in equations with state variables are the errors between the corresponding 
states of the reference system (system 1) and the process (system 2). We define the 
following set of variables 
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ei = Xi — X 2 , e 2 = Yi- Y 2 , e 3 = Z\ — Z 2 


The state equations of the plant in the error form are given by 


e x = -e 2 - e 3 - u 

(5-17) 

<22 = ei + ae 2 

(5.18) 

e 3 = e 3 (Zi - c) + ei(Zi - e 3 ) 

(5.19) 

This dynamical system representing the nonlinear plant has an c 

iquilibrium point at 

(0,0). The dynamics of the exo-system are assumed to be governed by the following 
equation 

iv = aw 3 , a < 0 

(5.20) 

The plant output (in error form) and the reference output are given 

by (Equation (5.3)) 

h{x) = 6l 

(5.21) 

q(w) = -w 2 

(5.22) 


The control objective now is to track the time varying reference signal from the exo- 
system such that the error e„ = (e, - w 2 ) -> 0 as t -+ 00. The exo-system has been 
selected such that it has a stable zero solution and the objective is to apply suitable 
control action in order to stabilise the zero solution of the plant. This would then 
ensure that the process tracks the reference trajectory. 

The linearized form of the composite system comprising of the plant and the exo- 
system has one zero eigenvalue. Hence this system has a center manifold, which is 
he graph of the mapping e = „( w ). The functions »(„) are evaluated by employing 

* Ce T, “ anif0ld the0rem ' T1 ' eSe fUn ° ti0nS are obtei '«“ b y solving the partial 
. eren isd equations ( 5 . 5 ). Computing the functions rr( w ) involves solving the partial 

Ifferentia equations. Instead we approximate the solution assuming the solution to 


n(w) - ciw 2 + c 2 w 3 + C3U/ 4 + C4W 5 
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Setting the error of the composite system equal to zero, we obtain the first compo- 
nent of 7r(w) to be e\ = ni (w) = w 2 . The other two components of n(w) are obtained 
by solving the partial differential equations as explained in the previous section. The 
coefficients of the polynomial approximation are given by 

c 2 = 0 
c 4 = 0 
c 2 = 0 
c 4 = 0 

The feedback control law which solves the regulator problem is now given by 

u = a(x,tu) = c(w) 4- &i(ei — w 2 ) + k 2 {e 2 - 7r 2 (w)) + 
k 3 {e 3 - r 3 (w)) 

c{w) = —2 ww — 7 r 2 (w) — tt 3 (w) 

We integrated the system (5.11)-(5.16) using the above control law. Results of the 
simulations are shown in Figs. 5.3(a)-(e). Fig. 5.3(a) and 5.3(b) show the variation 
in the error components e x — X\ — X 2 , e 2 = Yi — Y 2 and e 3 = Z\ — Z 2 with time 
respectively. Evolution of the differences (e x — w 2 ), (e 2 — n 2 (w)) and (e 3 — tt 3 (w)) with 
time is shown in Fig. 5.3(c) and 5.3(d) respectively. It can be seen that the evolution 
of the error between the system states corresponds to the evolution of the exo-system, 
thus rendering the center manifold invariant. The evolution of the control variable u 
is shown in Fig. 5.3(e). 

In the previous chapter our objective was to design a controller for controlling the 
process output along a desired chaotic trajectory. As we saw in that chapter, tracking 
was successful only when the zero dynamics of the system were asymptotically stable. 
The controller design discussed in this chapter is based on a more stringent requirement. 
Here, It is desired to design the controller for tracking the reference output, but with the 
additional requirement of an asymptotically stable linearized form of the closed-loop 
system. This is reflected in the form of the control law, which consisted of two distinct 


(5.23) 

(5.24) 


-1 


e 2 = ir 2 (w) : c x = - 


c 3 


2aci 


e 3 = 7T 3 (w) : Cl = 


Xi-c 


r — £ 1 + 2 acj 
03 Xi-c 



Figure 5.3: Evolution of the errors and the manipulated variable with time, (a): 
e i = (*i ~ x * and e 2 = (*i - la); (b): e 3 = (Z x - Z 2 ); (c): e cmI = {e x - 7Tx(u/)) and 
e C m2 = (e 2 - tt 2 («;)); (d): e cm3 = (e 3 - 7r 3 (u>)); (e) u. a = -1, /c x = 12, fc 2 = 4 and 
ft 3 = — 1 in (5.24). 
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Figure 5.4: Steady-state bifurcation diagram depicting the dependence of the steady- 
states of (a) Y\ and (b) Y 2 on r. 

parts (Equations (5.7) and (5.23)). The first part of the control law ensured tracking of 
the reference signal, while the second part ensures stability of the the linearized form 
of the closed-loop system. 

5.4 An Autocatalytic Reaction in a CSTR 

We now consider the cubic autocatalysis reaction with decay[8] 

A + 2 B — y 3 B 
B -y C 

in an isothermal CSTR. The dimensionless mass balance equations for the species A 
and B respectively are given by 

Yi = r(l - Yi) - Y x Y 2 2 (5.25) 

y 2 = T (f 3 -Y 2 ) + Y 1 Yi-kY 2 (5.26) 

The overdot indicates differentiation with respect to time. The dimensionless param- 
eters are defined as follows: r is a dimensionless characteristic time of the system, j3 
is the dimensionless feed concentration of species B and k is a dimensionless reaction 
rate constant. In this study, we fix /? = 0.10 and k = 0.080. The dependence of the 
steady-states of Yi and Y 2 on r is shown in Fig. 5.4. 
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There are two regions of t where multiple steady-states occui. Howuvci, the legion 
of interest to us, is the region where operating the reactor at two values ol r (the control 
input), gives the same steady-state value of the controlled output. Oi oi \ 2 ). In this 
section our objective will be to design a nonlinear controller to stabilize the system 
at its desired steady-state. As mentioned in the first section of this chapter, several 
differential-geometry control techniques are of limited usefulness in such situations 
where the system can be non-minimum phase for the given choice of the controlled 
output. 

For ease of application of the method described in section 5.2, we rewrite the above 
equations in deviation form, by defining new variables X and Y as follows 

X = X.-Y U Y = Y,~Y 2 

Here, X s and Y s represent the desired steady states of the system. In the new co- 
ordinates, the model equations (5.25) and (5.26) are given by 

AT = -r(l + X -X s )-(X 3 -X)(Y,-Y) 2 (5.27) 

Y = - r (P + Y-Y,)- (X, - X){Y 3 - Yf + k{Y, - F ) (5.28) 

We again select the exo-system as 

w = aw 3 , a = -10.0 (5.29) 

It is desired to control the system (5.27)-(5.28) at its steady state A r = Y = 0. The 
parameter r is the control parameter. In the first case, we consider X to be the 
controlled output. Therefore, we have (Equation 5.3) 

h(x) = X, q(w) = - w 2 

Setting the error of the composite system equal to zero, we obtain the first component 
of tt(w) to be X = tti (w) = w 2 . The other component of tt (w) is obtained by solving 
the partial differential equations as explained in. the previous section. We choose 

Y = TT2{w) = CiW 2 + C 2 W 3 + C3W 4 + C4W 5 
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Figure 5.5: Response of the system states X and Y to a (a) +10% change and (b) 
—10% change in the input r. 

The coefficients in the above polynomial can be calculated to be 

. _ Y}(P-Y s +l-X,)-X a Y?+kY s 

Cl “ X,Y?-2XsY,(P-Y,+l-X.)+k{l-X,) 


-2c 1 Y,(0-Y a +l-Xs)-X,c‘>(P-Ys+l-X s )+Y?+ciY s 2 +2X,Y,c 1 +2X,Y,c2-kc 1 -2aci(l-Xs)+2a(P-Y s ) 

C 3 — XsY?+k(l-X s )-2X,Y,(P-Y,+l-X a ) 

and C 2 = 0, c 4 = 0 

The feedback control law for regulating the plant output is now given by 

r = q ( x , w ) — c(w) + ki(X — w 2 ) + k 2 (Y - tt 2 (w)) (5.30) 

-2aw 4 + (X s - w 2 ){Y s - tt 2 (w)) 2 ^ 

c(w) = a +»*-*.) (5 - 31) 

We study the controllability characteristics of this system for r = 0.038. The 
system exhibits input multiplicity for this value of r. We consider the set-point to be 
the middle steady-state in Fig. 5.4. This corresponds to the steady-state values of 
X s = 0.7576 and Y s = 0.1103. The response of the system to a ±10% change in the 
input r is shown in Fig 5.5. 

We compare the results obtained using our output regulation method with those 
obtained using the input-output linearization procedure of Kravaris and Chung[67]. 
This technique aims at synthesising a feedback controller such that the input-output 
map is linearized. The linearization is obtained by seeking a transformed input v = 



no 


v(x,u) such that the dependence of plant output y on v is linear, i.e., t-he input-output 
relationship can be represented by a linear differential equation ol the hum 


m 


2 > 


<h 

dP 


v 


(5.32) 


The conditions for the existence of such a map are discussed by Kt aval is and C. hung[67]. 
These developments are based on the techniques discussed in chapter 4. 

Such a map which linearizes the relationship between the plant output y and the 
transformed input v is obtained by (i) multiplying the (r + 1) equations (4.5) by 0j , j = 
0, 1, • • • , r and (ii) summing the (r+1) equations thus obtained. The transformed input 
v which linearizes the input-output map is given by 


v = fah(x) + 0iLth{x) + • • • + (3 r L}h{x) + 0 r L % L r f l h{x)u (5.33) 


The dynamics of the v — y system are determined once the constants fy, j = 
0, 1, - - - , r, in (5.32) are specified. The transformed input v can be taken to be from 
an external PI controller which forces the output y to track the desired response yd- 
Hence 


\{Vd -y) + ~ [ (yd- y) 

I T e JO J 


(5.34) 


Using (5.33) and (5.34), we can obtain the actual nonlinear feedback control law 



I<e 


Xvd ~ y) + £ fo(Pd - y)] - E5«o ^jLpi(x) 
P r L s L r f l h{x) 


(5.35) 


The input-output linearizing differential equation (5.32) above, is of order r. In 
most cases r < n, and this equation therefore, accounts only for part of the closed-loop 
dynamics of the process. Control of the system is successful only when the dynamics 
of the remaining (n — r) states representing the internal dynamics of the system are 
stable. As discussed in the previous chapter, stability of the internal dynamics is given 
by the stability of the zero dynamics. Consequently, this technique can be successfully 
used only for minimum-phase systems. 

In the case of the autocatalytic system, our objective is to regulate the system at 
its steady-state, i.e., at yd = 0. The system is of relative order r = 1, since we consider 
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Output 

Case 

K e 

r e 

fio 

X 

+20% 

5.0 

0.25 

-100.0 

X 

-20% 

0.80 

1.60 

-2.0 

Y 

+20% 

5.00 

0.50 

-10.0 

Y 

-20% 

10.0 

2.50 

-2.50 


Table 5.1: Parameter values. Pi = 1, for all cases. 

r to be the manipulated input. For this case, the linearized input-output map is of the 
form 

A * =t ' +0oV (5 - 36) 

The parameters P 0 and Pi are selected such that the eigenvalues of the system (5.36) 
are in the left-half plane. Specifying these constants and the parameters K e and r e 
in (5.34) completely determines the dynamics of the linear v — y system. The actual 
nonlinear feedback control law is obtained by solving (5.36) and (5.27) or (5.28) for r . 
We have used this technique to track changes in the set-point of the system. Table 5.1 
gives the values of the parameters used in the simulation studies. 

The output regulation technique can be considered to be a general form of the 
input-output linearization procedure. This is because in the latter, the external system 
(Equation (5.2)) and the reference output q(w) (Equation (5.3)) are of the form 

Piw = Pow + v 
q(w ) = w 

As mentioned earlier, Kravaris and Chung consider v to be an external input controlled 
by a PI controller. This input can be chosen to be be equal to zero. Comparing the 
above equations with (5.29) and (5.3) respectively, it can be seen that defining the 
external system to be linear, results in a linear relationship between the input and 
the output. A disadvantage of this scheme is that it is useful only for the control 
of minimum-phase systems. Selecting the external system to be nonlinear and such 
that its linearization has all eigenvalues on the imaginary axis, enables us to employ 
the center manifold theorem in the design of the controller. Application of this result 
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Plant output 

Case 

Magnitude 

k\ 

h 

X 

SP 

+20% 

1.0 

10.80 

X 

SP 

-20% 

1.0 

6.0 

X 

UMD 

+15% in P 

0.10 

6.0 

X 

UMD 

-5% in p 

1.0 

10.0 

X 

MD 

+15% in p 

0.40 

5.0 

Y 

SP 

+20% 

1.0 

0.20 

Y 

SP 

-20% 

4.0 

0.80 

Y 

UMD 

+20% in P 

5.0 

0.20 

Y 

UMD 

-5% in p 

6.0 

0.50 

Y 

MD 

+20% in P 

8.0 

0.20 

Y 

MD 

-5% in p 

6.0 

0.50 


lable 5.2: Parameter values used in the output regulation technique. SP: Change in 
set-point; UMD: Unmeasured disturbance; MD: Measured disturbance. 

enables us to succesfully apply the control technique for the trajectory tracking of 
non-minimum phase systems. 

The response of the system using the input-output linearisation procedure of Kravaris 
an Chung[67], for a ±20% change in the setpoint, is shown in Fig. 5.G. The system 

* S n0n - mmlmum the choice of the controlled output. It can be seen from the 

gures that the system has stabilized at the steady-state corresponding to a non-zero 
ue o . In both cases, the manipulated input has stabilized at a value close to 

reaulaf 16 T" ** * ±2 ° % “ “» “‘point using the output 

in the r n ^ T in Kg ' 5 ' 7 ' ThC Va ‘ UeS ° f thc para,Mtos *« ««I h used 
to be 2 Z m TablC 5a ' Applyh,g this «*** can cause the input r 

from th f TC ' T T 0 ^ MS ^ reSettin8 ' " °' 38 - whenever , < 0. It is seen 
The rest) 168 1™ reSettmg affects the evolutio11 of X more than it alfects Y. 
shown pTTT °H t SyStem U " meaSUred “ d disturbances in fi is 

of the relative order oUhelyZl^th' ^ 4 ’. ”* PreSented the defi “ ti0 ” 

respect to the input. The relative order of 
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Time 


Figure 5.6: Response of the system states ((a), (c)) and the input ((b), (d)) to a +20% 
change ((a), (b)) and —20% change ((c), (d)) in the set-point using the input-output 
linearization technique. 
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Figure 5.7: Response of the system states ((a), (c)) and the input ((b), (d)) to a 

+20% change ((a), (b)) and -20% change ((c), (d)) in the set-point using the output 
regulation technique, w = w 2 . 
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Figure 5.8: Response of the system states ((a), (c)) and the input ((b), (d)) to unmea- 
sured disturbances of +15% ((a), (b)) and —5% ((c), (d)) in /?. w = w 2 . 

the system with respect to the disturbance can be defined along similar lines. For the 
case under consideration, the system is of relative order one with respect to the input 
and of relative order two with respect to the disturbance. Somewhat surprisingly, it 
is observed that the disturbance affects X more than it affects Y. The disturbance 
results in an off-set in the steady-state response of the system. It was observed that 
the off-set increased rapidly on increasing the magnitude of the disturbance. 

We next consider the case when Y is the controlled output. Therefore, we have 
(Equation 5.3) 

fi(x) = Y, q(w ) = — w 2 

Setting the error of the composite system equal to zero, we obtain the first component 
of 7 t(w) to be y = ni(w) = w 2 . The other component of tt(w) is obtained by solving 
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Figure 5.9: Response of the system states ((a), (c)) and the input ((b), (d)) to measured 
disturbances of +15% ((a), (b)) and -5% ((c), (d)) in (3. w = w 2 . 
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the partial differential equations as explained in the previous section. We choose 

X = 7 r 2 (w) = CiW 2 + c 2 w 3 + C 3 W 4 + C 4 W 5 

The coefficients in the above polynomial can be calculated to be 

„ - -2x M Y,(i-x,))+k(i-x t )+x.Y,*Y'-2X,Y M (p-Y.) 

1 - Y*(l-X.yx,Y*+kY,+Y*(p-Y„) 


-2c 1 a(P-Ys)+2a(l-X,)+kci-Y,Y s ^+2Y.ci*(l~X,)+Xs(.l-X.) 

Y?(l-X.)-X,Y?+kY,+Y?tf-Ys) 


2X,Y,c 1 +2Y s c 1 (P-Y,)-Y}c 1 +XAl3-Y 3 )-2X,Y s 

Y?(l-X t )-X.Y?+kY»+Y*(P-Y.) 


and c 2 = 0, e 4 = 0 

The feedback control law for regulating the plant output is now given by 

r = oc(x,w) = c(w) + (ki(Y — w 2 )) + (k 2 (X — ir 2 (w))) (5.37) 

-2aw 4 - {X s -tt 2 (w)){Y s -w 2 ) 2 + k(Y s -w 2 ) 

C{W) = U>-Y. + rf) (5 - 38) 

We again study the controllability characteristics of this system for r = 0.038. The 
response of the system using the input-output linearization procedure of Kravaris and 
Chung[67], for a ±20% change in the setpoint, is shown in Fig. 5.10. The results 
indicate that the system can be non-minimum phase (for the +20% change in the set- 
point) for this choice of the controlled output. The response of the system for a ±20% 
change in the set-point using the output regulation technique is shown in Fig. 5.11. 
The values of the parameters k\ and k 2 used in the simulations are listed in Table 5.1. 

The response of the system to unmeasured and measured disturbances in-/? is shown 
Fig. 5.12 and Fig. 5.13 respectively. For this case of the controlled output of the plant 
(y), the system is of relative order one both with respect to the input as well as 
the disturbance. The disturbance results in an off-set in the steady-state response of 
the system. The results are similar to those obtained with the case when X was the 
controlled output. 
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Figure 5.10: Response of the system states ((a), (c)) and the input ((b), (d)) to a +20% 
change ((a), (b)) and -20% change ((c), (d)) in the set-point using the input-output 
linearization technique. 
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Figure 5.11: Response of the system states ((a), (c)) and the input ((b), (d)) to a 
+20% change ((a), (b)) and -20% change ((c), (d)) in the set-point using the output 
regulation technique, w = w 2 . 
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Figure 5.12: Response of the system states ((a), (c)) and the input ((b), (d)) to un- 
measured disturbances of +15% ((a), (b)) and -5% ((c), (d)) in /?. w = w 2 . 
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Figure 5.13: Response of the system states ((a), (c)) and the input ((b), (d)) to mea- 
sured disturbances of +15% ((a), (b)) and —5% ((c), (d)) in /?. ui ~ w 2 . 
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5.5 Conclusions 

It is well known that not all chaotic systems can be synchronized with the presently 
available methods of synchronization. The Rossler system is an example of such a 
system. We have employed recent results from nonlinear control theory and presented 
a method for the synchronization of such systems. 

The method has been applied to develop a control technique for an isothermal 
reactor system exhibiting input multiplicity. The results obtained using the output 
regulation technique have been compared with those obtained using the input-output 
linearization procedure. 

A limitation of the method appears to be its inability to handle models containing 
exponential nonlinearities. This is because applying the center manifold theorem to 
models containing these nonlinearities requires that we approximate the exponential 
term with the first few terms of its Taylor series expansion. This can lead to a drastic 
reduction in the performance of the controller. 
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Appendix A 

Singularity Theory 

The steady states of lumped-parameter nonlinear systems are determined by solving 
systems of nonlinear algebraic equations. These systems are often characterized by 
a large number of parameters. In such cases it can be very difficult to determine the 
behavior of the system by an exhaustive study. In many situations of practical interest, 
two important issues which need to be addressed are 

• What is the number of possible steady-states the system can possess and 

• How are the steady-states influenced by changes in one or more of the system 
parameters? 

Application of singularity theory enables us to answer these questions in an elegant 
and comprehensive manner. 

In this Appendix, we summarise the features of singularity theory [35], with empha- 
sis on the applicability of the technique to the coupled reactor-separator system. 

We begin by assuming that at steady-state, the model of the system under consid- 
eration can be represented by a single nonlinear algebraic equation in one of the state 
variables x 

f(x,\,p)=0 (A.1) 

Here x represents the dependent (state) variable, A represents the control parameter 
and p is the parameter set containing all other parameters. We are interested in 
studying the dependence of x on A. 
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The multiplicity features of the system depend on the choice of the parameters in 
the parameter set p. These are represented as bifurcation diagrams. A bifurcation 
diagram describes the steady-state dependence of the state variable x on the control 
parameter A. This control parameter is known as the bifurcation parameter. The 
bifurcation diagrams are plotted by fixing all the other parameters p. In this discussion, 
a bifurcation point is defined as a point in parameter space, where there is a change in 
the number of solutions. 

The multiplicity features can change only when critical surfaces in the parameter 
space p are crossed. Singularity theory enables us to identify these surfaces and to 
classify the parameter space into regions having different multiplicity features. The 
bifurcation diagrams in each region of the parameter space are therefore qualitatively 
different. 

A detailed presentation of the mathematical theory is given in Golubitsky and Scha- 
effer [35]. Balakotiah and Luss[22-24] have applied the theory to analyse the steady- 
state multiplicity features of several lumped-parameter systems. 

The first step in the application of singularity theory consists of identifying the 
normal form of the system. This is the simplest polynomial g(y, p) = 0, depending 
on a parameter p, to which the system /(x, A) is equivalent. The notion of equivalence 
is defined by conditions for the two solution sets of the equations / = 0 and g — 0 
to be qualitatively similar. These conditions reduce to computing the values of the 
some of the parameters from the parameter set p, at which the system / = 0 has 
an identical qualitative behavior as the polynomial g = 0 at the origin. Singularity 
theory provides a list of the calculations which must be performed to recognize a given 
equation /(x, A) to be of a given qualitative type. The point in parameter space at 
which the system is equivalent to the simple polynomial g(y, p) is the singular point 
of the system. If the function / satisfies the conditions of equivalence at the singular 
point in parameter space p, then we expect the solution sets of the system / = 0 to 
be similar to the solution sets of the polynomial g — 0 in the neighborhood of this 
point. This follows since / will have the same Taylor series expansion at the singular 
point as the polynomial g has about the origin. A list of normal forms is tabulated in 
Golubitsky and Schaeffer [3 5]. 
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Figure A.l; Bifurcation diagrams showing the change in the multiplicity features (a) 
before crossing the hysteresis variety; (b) on the hysteresis variety and (c) after crossing 
the hysteresis variety. 


parameters in p are kept constant. The hysteresis variety can thus be represented 
as a curve in a two-parameter plot. For parameter values on the hysteresis variety, 
the bifurcation diagram shows a point of inflexion or a point with a vertical 
tangent. Across the variety, two new solutions appear or disappear (depending 
on the direction of crossing). Fig. A.l shows a schematic of the bifurcation 
diagram on the hysteresis variety and on either side of it. 


Isola Variety: This critical surface is generated by solving 


o 

1 dx d\ 


(A.3) 


These equations are solved for x, A and a parameter from the parameter set p, 
while keeping the others constant. Crossing the isola variety can either result in 
an isola bifurcation or a simple bifurcation. In the case of an isola bifurcation, 
an isolated branch of solutions either appear or disappear as the isola variety is 
ciosscd. In the case of a simple bifurcation, two solution branches merge and 
then separate. These are shown in Fig. A.2. 


When there are no constraints acting on the state variable x, multiplicity features of the 
system can change only when the hysteresis or the isola varieties are crossed. However, 
as we have seen in the third chapter, steady-state behavior of the reactor-separator 
system can be influenced by the presence of constraints on the conversion z. To recall, 
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Figure A. 2: Bifurcation diagrams showing the change in the multiplicity features (a) 
before crossing the isola variety; (b) on the isola variety and (c) after crossing the isola 
variety. Figures in the first and second rows are for the isola bifurcation and the simple 
bifurcation respectively. 

feasible steady-state operation of the system is possible only when y e < z < x e . Such 
constraints on the state variable appear as feasibility boundaries in the bifurcation 
diagrams. New types of multiplicity features arise depending on the manner in which 
these boundaries are crossed. Singularity theory can also be employed to account 
for these new kinds of multiplicity features in a systematic way. Special parameter 
sets have been defined to account for the presence of the feasibility boundaries[24, 
35]. These enable us to demarcate the locus of the feasibility boundaries in a suitable 
parameter space. We now summarise these special parameter sets. We have retained 
the notation used in section 3.3 of the third chapter. 

• Boundary limit set (BLS): The parameter set corresponding to which a turn- 
ing point of the steady-state branch in the bifurcation diagram lies on the fea- 
sibility boundary forms the BLS. This set is given by the solution to the set of 
equations 

G = G z = 0 at z = x e or z = y e (A.4) 
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• Boundary tangent set (BTS): For parameter values belonging to this set, 
the steady-state branch in the bifurcation diagram is tangential to the feasibility 
boundary. This is obtained by solving 

G = Gr = 0 at z = x e or z = y e (A.5) 

• Double cross set (DCS): This is given by the values of x e - y e for which two 
branches of the bifurcation diagram cross the feasibility boundaries simultane- 
ously. This is obtained by the solution to the set of equations 

G(z = x e ) = G(z = y e ) — 0 (A.6) 

In chapter 3, these equations have been solved for r and (3q for fixed values of all the 
other parameters. Fig. A.3 shows the different types of bifurcation diagrams resulting 
due to the presence of these feasibility boundaries. 

In section 3.4 of the third chapter, we have considered ^ to be the bifurcation 
parameter for the case of isothermal-isobaric flash and when F and F 0 are the specified 
variables. In this case, the feasibility boundaries (BLS) appear on the bifurcation 
parameter 2 . Schematic representation of the bifurcation diagrams are shown in Fig. 
A.4. 

Multiplicity features for the case of isothermal and isobaric separation when (Mr, F) 
are the specified variables (section 3.3.2), can change when another critical surface is 
crossed in the parameter space. This is the limit and cross set (LCS). Fig. A.5 shows 
the different types of bifurcation diagrams when this set is crossed. This set is obtained 
by solving the following sets of equations 

G(z = 3 /e) = 0, G = G z = 0 (A. 7) 

In section 3.1 of chapter 3, these equations have been solved for z , r and (3q for fixed 
values of all the other parameters. 

In section 3.3 of the third chapter, we have applied singularity theory to determine 

• all the possible local bifurcation diagrams and 

• classifying a suitable parameter space into regions having different multiplicity 
features 

for different sets of specifications of the reactor-separator system. 
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Figure A.3: Bifurcation diagrams showing the change in the multiplicity features in 
the presence of feasibility boundaries, (a) before crossing the boundary sets; (b) on 
the boundary sets and (c) after crossing the boundary sets. Figures in the first, second 
and third rows are for the BLS, the BTS and the DCS respectively. 







Figure A.4: Bifurcation diagrams showing the change in the multiplicity features for 
the case of (7>, P) and (F, F 0 ) flash, (a) before crossing the BLS; (b) on the BLS and 
(c) after crossing the BLS. 
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Figure A. 5: Bifurcation diagrams showing the change in the multiplicity features for 
the case of ( Tp,P ) and ( Mr,F ) flash, (a) before crossing the LCS; (b) on the LCS 
and (c) after crossing the LCS. 




